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We study the zero temperature superfluid-insulator transition for a two-dimensional model of 
interacting, lattice bosons in the presence of quenched disorder and particle-hole symmetry We 
follow the approach of a recent series of papers by Alt man, Kafri, Polkovnikov, and Refael, in which 
the strong disorder renormalization group is used to study disordered bosons in one dimension. 
Adapting this method to two dimensions, we study several different species of disorder and uncover 
universal features of the superfluid-insulator transition. In particular, we locate an unstable finite 
disorder fixed point that governs the transition between the superfluid and a gapless, glassy insulator. 
We present numerical evidence that this glassy phase is the incompressible Mott glass and that the 
transition from this phase to the superfluid is driven by a percolation- type process. Finally, we 
provide estimates of the critical exponents governing this transition. 



I. INTRODUCTION 

In the 1980s, seminal experiments on helium adsorbed 
in Vycor first attracted the attention of theorists to the 
random boson problem 1 ' 2 . The onset of superfluidity in 
this disordered system showed, in some respects, strik- 
ing similarities to ideal Bose gas behavior. Thus, in the 
decade before the realization of Bose-Einstein conden- 
sation in cold atomic gases, disordered bosonic systems 
were actually proposed as possible realizations of this elu- 
sive phenomenon. While studies of disordered bosons did 
not ultimately lead to the observation of Bose-Einstein 
condensation, the random boson problem continued to 
stimulate theoretical activity because of its considerable 
richness. The interplay of interactions and disorder leads 
to a variety of phases in bosonic systems; the description 
of these phases and the transitions between them is an 
ongoing challenge 3 ' 4 . 

The theoretical investigation of random bosons in one 
dimension was pioneered by Giamarchi and Schulz, who 
described the transition to superfluidity in the presence 
of perturbatively weak disorder 5 . Subsequently, Fisher, 
Weichman, Grinstein, and Fisher expanded upon the 
work of Giamarchi and Schulz by establishing the now 
canonical zero temperature phase diagram of the Bose- 
Hubbard model with chemical potential disorder. The 
superfluid and Mott insulating phases of the clean model 
are separated by another insulating phase, the Bose glass. 
In this glassy phase, there exist rare-regions in which the 
energetic gap for adding another boson to the system 
vanishes. However, any additional bosons are localized 
by the disordered environment, rendering the phase a 
gapless, compressible insulator. Fisher et al. argued that 
there should be no direct superfluid-Mott insulator tran- 
sition in the presence of quenched, uncorrelated disorder, 
or in other words, that the Bose glass always intervenes 
between the phases of the clean model 4 . The general pic- 
ture formulated by these authors has been vindicated by 
over two decades of subsequent theoretical and numerical 
work 6 . On the experimental front, following the observa- 
tion of the Mott insulator to superfluid transition in cold 



atomic gases by Greiner et al. 7 , there has been progress 
towards the realization of a Bose glass, including sugges- 
tive but inconclusive evidence that the phase has already 
been observed 8 11 . 

Meanwhile, evidence has accumulated that disordered 
bosonic systems may exhibit more exotic glassy phases in 
the presence of particle- hole symmetry. One such phase, 
the so-called Mott glass, differs from the Bose glass in 
that it is incompressible. This phase was originally pro- 
posed for systems of disordered fermions by Giamarchi, 
Le Doussal, and Orignac, but these authors predicted 
that it can exist in bosonic systems as well 12 . Subse- 
quently, Altman, Kafri, Polkovnikov, and Refael studied 
a variant of the Bose-Hubbard model in which chemical 
potential disorder is omitted in favor of strong disorder 
in the on-site interaction and hopping. In the limit of 
large, commensurate boson filling, this model is equiv- 
alent to a chain of quantum rotors that can describe a 
Josephson junction array. Relative to the large filling, 
this rotor model exhibits an exact particle-hole symmetry 
which has important consequences for the phase diagram. 
Through a strong disorder renormalization group analy- 
sis, Altman et al. found that this symmetry results in the 
appearance of a Mott glass phase in this model 13 . The 
same authors later considered a generalization of this ro- 
tor model which allows for random offsets in the filling, 
effectively reintroducing a chemical potential. Generic 
disorder in the random offsets violates the particle-hole 
symmetry, and in this case, the rotor model exhibits the 
usual Bose glass phase 14 ' 15 . This breakdown of the Mott 
glass demonstrates the link between exotic phases and 
symmetries of the disordered Hamiltonian. Irrespective 
of the identity of the glassy phase, Altman et al. also 
found that the superfluid-insulator transition at strong 
disorder lives in a different universality class from the 
weak disorder transition of Giamarchi and Schulz 13 ' 14 . 

In this paper, we study the two-dimensional analogue 
of the rotor model considered by Altman et al. in order 
to investigate the particle-hole symmetric random boson 
problem in d > 1. One of our goals is to look for the 
Mott glass, but more generally, we aim to characterize 
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FIG. 1: A schematization of the universal features of the pro- 
posed flow diagram. The x-axis gives the ratio of the mean 
of the renormalized Josephson coupling distribution to the 
mean of the renormalized charging energy distribution. The 
y-axis gives the ratio of the standard deviation of the Joseph- 
son coupling distribution to its mean. In this context, the 
Josephson coupling distribution only includes the dominant 
2N couplings in the renormalized J distribution, where N is 
the number of sites remaining in an effective, renormalized 
lattice. See the text of Section III for the reasoning behind 
the exclusion of weaker Josephson couplings from statistics. 



the phases and the superfluid-insulator transition of our 
model. Like Altman et al., the tool we use to do this 
is the strong disorder real space renormalization group 
(RG), first formulated by Ma and Dasgupta 16 to study 
the ID Heisenberg antiferromagnet about thirty years 
ago. A numerical application of the RG by Bhatt and 
Lee followed shortly thereafter 17 , and the method was 
later expanded upon and applied to more general spin 
models by Daniel Fisher 18 . Strong disorder renormaliza- 
tion has proven to be a powerful tool in the analysis of 
several disordered systems, especially in one dimension. 
In higher dimensions, application of the RG has been 
rarer, because in addition to the generic intractability 
of analytical approaches 19 , there are few known transi- 
tions that exhibit so-called infinite randomness, a prop- 
erty that guarantees that the RG becomes asymptoti- 
cally exact near criticality. The random boson model is 
not expected to exhibit infinite randomness, and indeed, 
the numerical data we present in this paper is consistent 
with this expectation. Hence, in addition to physical 
questions about the phases and phase transitions of the 
model, our work also aims to address a methodological 
question: might the strong disorder RG give useful infor- 
mation about a model, even when confronted with the 
twin difficulties of higher dimensionality and the absence 
of infinite randomness? 



Summary of the Results 

Our main results are as follows: we present numerical 
evidence for the existence of an unstable finite disorder 
fixed point of the RG flow, near which the distributions 
of Josephson couplings Jjk and charging energies Uj in 
the rotor model flow to universal forms. A schematic 
picture of this unstable fixed point and the flows in its 
vicinity is given in Figure 1. 

To the left of the diagram, flows propagate towards a 
regime in which the ratio of J, the mean of the Joseph- 
son couplings, to J7, the mean of the charging energies, 
vanishes; meanwhile, the ratio of AJ, the width of the 
Josephson coupling distribution, to J grows very large. 
These flows terminate in one of two insulating phases. 
The first is a conventional Mott insulator, in which it 
is energetically unfavorable for the particle number to 
fluctuate from the large filling at any site. The other is 
a glassy phase, in which there exist rare-regions of su- 
perfluid ordering. As the thermodynamic limit is ap- 
proached, arbitrarily large rare-regions appear, driving 
the gap for charging the system to zero. However, the 
density of the largest clusters decays exponentially in 
their size, and the size of the largest cluster in a typi- 
cal sample does not scale extensively in the size of the 
system. Moreover, the largest clusters are so rare that 
they cannot generate a finite compressibility. Thus, the 
phase is a Mott glass. 

This insulating phase gives way to global superfluidity 
when the rare-regions of superfluid ordering percolate, 
producing a macroscopic cluster of superfluid ordering. 
The appearance of the macroscopic cluster is associated 
with flows that propagate towards the lower right of Fig- 
ure 1, indicating that the unstable fixed point governs 
the glass-superfluid transition. Our numerical implemen- 
tation of the strong disorder RG allows us to extract es- 
timates for the critical exponents that characterize this 
transition. We are thus able to construct a compelling 
picture of the superfluid-insulator transition: a picture 
that must, however, be checked by other methods be- 
cause of the perils of employing the strong disorder RG 
method in the vicinity of a finite disorder fixed point. 



Organization of the Paper 

We begin our analysis of the disordered rotor model 
in Section II by introducing the model, noting its rela- 
tionship to the standard disordered Bose-Hubbard model 
and its special symmetries. We also discuss the clean 
limit and the disordered problem in one dimension. Sec- 
tion III is devoted to a description of the strong disor- 
der renormalization group, as applied to the disordered 
rotor model. In discussing the method, we emphasize 
the adaptations needed to use the technique in dimen- 
sions greater than one. We then present data collected 
from our numerical implementation of the strong disor- 
der renormalization procedure in Section IV and sub- 
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sequently explore what the data can tell us about the 
zero temperature phases and quantum phase transitions 
of our random boson model in Section V. In Section VI, 
we summarize the results, make connections to experi- 
ments, and give an outlook. 

II. THE MODEL 



The operators fij now correspond to the particle num- 
ber deviations from the large filling No. As such, nj can 
take on any integer value from — Nq to oo, but we assume 
that No is so large that we can let nj run from — oo to 
oo. The same approximation allows us to drop sublead- 
ing (in -^-) terms in the hopping. We finally define the 
couplings Jjk = 2tjkNo to arrive at the quantum rotor 
Hamiltonian: 



In motivating the model that we study in this paper, 
we begin by writing down a disordered Bose-Hubbard 
Hamiltonian that includes randomness in the interaction 
and hopping along with the usual chemical potential dis- 
order: 

H hh = tjk(a]a k + a\aj) + ^ Ujd)dj(d)dj - 1) 

Ok) 3 

-Y^H^jCLj (1) 

3 

Here, the creation and annihilation operators satisfy 
bosonic commutation relations: 

[aj,a{] = 5 jk (2) 

and the hopping is between all nearest-neighbor sites on 
an L x L square lattice with periodic boundary condi- 
tions. An alternative representation of this model is given 
by the number and phase operators: 

aj = yJWj (3) 



H = - ^ J jk C0S(<^ - 4> k ) + ^ U 3^ 2 j ( 7 ) 
(jk) 3 

This model, constructed as the large filling limit of a 
Bose-Hubbard Hamiltonian, can also be viewed as a de- 
scription of a two-dimensional array of superconducting 
islands connected by Josephson junctions 13 ' 14 . Moreover, 
Vosk and Altman have recently demonstrated that the 
one-dimensional model is relevant to cold atomic gases 
of rubidium-87 20 . 

When the Josephson couplings Jj k and charging en- 
ergies Uj are uniform, the rotor model (7) exhibits a 
quantum phase transition between superfluid and Mott 
insulating phases at zero temperature. This transition is 
in the universality class of the three-dimensional classical 
XY model 3 ' 4 ' 21 , and one recent study determines that the 
transition occurs at jj ~ 0.345 22 . The critical exponent 
governing the divergence of the correlation length at the 
clean transition is v ~ 0.663 23 . This exponent violates 
the Harris criterion: 

vd > 2 (8) 



[4>j,N k ] = i&jk 
In terms of these operators: 



^bh — — ^2 

(3k) 



(4) 



+ ^N k e~ l(pk e l(p i^N j 



(5) 



To obtain a large, commensurate boson filling 7V , the 
chemical potential is tuned such that the on-site inter- 
action and chemical potential terms are minimized for 
Nj = Nq. This consideration fixes jij = (2 No — l)Uj. 
Then, if we expand the number operators around this 
large filling as Nj = Nq + n^, the Hamiltonian becomes: 



uk) 



i + 



nk_ 
N 



+ A 1 



n k 



No 

+ E Ujftj + (const.) 



1 + 



No 



(6) 



when d = 2. Violation of the Harris criterion generically 
indicates that disorder will either change the universal- 
ity class of the clean transition or completely smear it 
away. In the former case, one or more Griffiths phases 
will separate the phases of the clean system 24 ' 25 . The 
nature of this intervening region depends upon the spe- 
cific model in question. In ID at T = 0, Altman et al. 
found an incompressible Mott glass phase and a quantum 
phase transition of Kosterlitz-Thouless type between this 
glassy phase and the superfluid 13 . The RG fixed point 
that controls this transition actually occurs at a point 
in the flow diagram where all Uj =0, meaning that the 
transition can be tuned by only varying the disorder in 
the Josephson couplings at arbitrarily weak interaction 
strength. 

In our work, we introduce disorder by choosing the 
initial distributions of charging energies and Josephson 
couplings, Pi(U) and Pi(J) to have one of the following 
forms: 

1. Gaussian distributions truncated at three standard de- 
viations: 



Pi(x) oc exp 



(x - x ) 2 



2al 



(9) 



for x G (xo — 3<7 X , xq + 3cr x ). 
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2. Power law distributions with upper and lower cutoffs: 



Pi{x) 



T] + 1 



77+1 _ »+r 



(10) 



for X G (^min) ^max)* 

3. Flat distributions with upper and lower cutoffs: 
= 



(11) 



u max ^min 



for x G (#min?#max)- Of course, this is just a power 
law with exponent r] = 0. 

4. "Bimodal" distributions consisting of two flat peaks 
centered at xg and x^: 



Pi(x) = 



2Sx 



(12) 



for x G (x£ — ^r,xi + 4f ) and x G (x^ — ^f,x h + 4f ). 

All of these distributions have positive lower and upper 
cutoffs (x m i n and x max respectively) and have zero weight 
for x outside of these bounds. This restriction avoids 
the complications of frustration in the phase degrees of 
freedom and the pathologies of the particle sinks that re- 
sult from on-site charging spectrums that are unbounded 
from below. 

Even in the presence of disorder, the Hamiltonian (7) 
respects two important symmetries. First, there is the 
global U(l) phase rotation symmetry: 



(13) 



This means that the Hamiltonian conserves total particle 
number: 



rhot 



(14) 



The model is also globally particle-hole symmetric: 



(15) 

The particle-hole symmetry exists because the chemical 
potential coupling to the true particle number Nj has 
been tuned precisely to the value that enforces the large, 
commensurate filling. If this chemical potential is al- 
lowed to deviate from this value, then it would manifest 
in a chemical potential coupling to the particle number 
fluctuation rij, or equivalently, in offsets to the large fill- 
ing. Such terms are absent from our rotor Hamiltonian 
(7), but let us momentarily consider the on-site charging 
spectrum (as a function of rij) in the general situation 
where a filling offset Srij may be present: 



The integer value of rij that minimizes this energy 
changes at half-integer Srij . For sufficiently strong disor- 
der in Srij , the introduction of an arbitrarily small global 
chemical potential shift will bring a finite fraction of sites 
j arbitrarily close to these density changing points. Thus, 
a finite density of particles will be added to the system, 
making it compressible. Nevertheless, these particles will 
be localized by the disordered environment, leaving the 
system globally insulating. This is the mechanism behind 
the formation of the Bose glass 4 . The situation changes 
at two special particle-hole symmetric points. One oc- 
curs when Srij is restricted to be integer or half-integer. 
Then, at the half-integer sites, there is a degeneracy in 
the charging spectrum: 



Srij 



1 



(17) 



E j( n j) = u Mi - 8n i f 



(16) 



In the ID model, Altman et al. showed that this degen- 
eracy gives rise to a random singlet glass 14 , but we will 
not explore this situation in this paper. Instead, we will 
focus on the other particle-hole symmetric point where 
Srij = for all sites j. Now, if the Uj are distributed such 
that U m i n > 0, the on-site charging spectrum always has 
a unique minimum that is protected by a gap. The usual 
mechanism for Bose glass formation is evaded, and this 
allows for the possibility of realizing more exotic glassy 
phases. Thus, the particle-hole symmetry is a crucial 
feature of our model that can influence the nature of the 
intervening glassy phase. 



III. STRONG DISORDER RENORMALIZATION 
OF THE 2D DISORDERED ROTOR MODEL 

As mentioned previously, the tool that we use to study 
the disordered rotor model is the strong disorder real 
space renormalization method. Here we briefly review 
the method and discuss its application to the model at 
hand. 

At first glance, problems involving strong quenched 
disorder may appear to be substantially more compli- 
cated than their clean counterparts. However, one way to 
motivate the strong disorder renormalization procedure 
is to consider that, in some cases, strong disorder can ac- 
tually serve as an advantage. In particular, disorder can 
make the problem of finding the quantum ground state of 
a model more local. Having identified the strongest of all 
the disordered couplings in the Hamiltonian, we can then 
use the assumption of strong disorder to argue that other 
couplings in the vicinity (in real space) of this dominant 
coupling are likely to be much weaker. This means that 
the ground state can locally be approximated by satisfy- 
ing the dominant coupling. Other terms in the Hamil- 
tonian can then be incorporated as corrections. Quite 
often, these other terms are treated by means of per- 
turbation theory, but this need not always be the case. 
These corrections manifest as new couplings in the model, 
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and thus, the procedure yields a new effective Hamilto- 
nian. Since part of the ground state is specified in this 
step, some degrees of freedom of the system are deci- 
mated away. By repeating the procedure, we can itera- 
tively specify the entire ground state 16 ' 17 . Moreover, we 
can examine the way in which the probability distribu- 
tions of the disordered couplings flow as the renormal- 
ization proceeds. One possibility is that the the model 
looks more and more disordered at larger length scales 
near criticality. This is the case for random transverse 
field Ising models in one and two dimensions 19,26,27 . In 
such cases, the model flows towards infinite randomness, 
and the strong disorder renormalization group becomes 
asymptotically exact near criticality 18 . Of course, it is 
also possible to flow towards finite or weak disorder, and 
in these cases, the reliability of the RG is less clear. We 
discuss this issue extensively, as it pertains to the disor- 
dered rotor model, in Appendix D. 

A. The Basic RG Steps 

We now concretize these ideas by application to the 
rotor Hamiltonian (7). In our model of random bosons in 
2D, there are two types of disordered couplings, namely 
charging energies and Josephson couplings. In each step 
of the renormalization, we identify the maximum of all 
of these couplings, which defines the RG scale: 

Sl = max[{U j },{J jk }] (18) 

How we then proceed depends upon which type of cou- 
pling is dominant. 




FIG. 2: The site decimation RG step: The site marked with 
the X has the dominant charging energy and is decimated 
away, generating bonds between neighboring sites j and k 
with the effective coupling given in equation (21). The new 
local structure of the lattice is shown to the right. 
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In the terms giving the perturbative corrections, we as- 
sume that the number fluctuations on all neighboring 
sites except m remain unmodified from their values in 
{n k }. We next consider the matrix elements of these 
states in Hx- Up to a constant term, these matrix el- 
ements are identically those that would result from an 
effective Josephson coupling: 

j JjxJxk f9U 
Jjk - — ^- — (21) 

between each two sites that were coupled to site X before 
the decimation step 13 . This process of site decimation is 
illustrated in Figure 2. 



Site Decimation 

Consider the case where the charging energy on site 
X is dominant. We define a local Hamiltonian in which 
this charging energy term is chosen to be the unperturbed 
piece. All Josephson couplings entering the correspond- 
ing site are considered to be perturbations: 

H x = U x n x - ^2 Jxk cos ^ x ~ ^ ( 19 ) 

k 

Satisfying the dominant coupling means setting tlx = 
to leading order. This defines a degenerate manifold of lo- 
cal ground states: |0, {n k }). In these kets, the first term 
corresponds to zero number fluctuation on site X and 
the second specifies the number fluctuations on all sites 
connected to X by a Josephson coupling. The degener- 
ate space is infinitely large, corresponding to all possible 
choices of {n k }. However, all matrix elements of the per- 
turbative Josephson couplings in this ground state mani- 
fold vanish. The leading corrections then come from sec- 
ond order degenerate perturbation theory, in which we 
calculate corrections coming from excited states: 

\0,{n k })' « |0,{n fe }) (20) 



Link Decimation 

Now, suppose that a Josephson coupling is the dom- 
inant energy scale. In this case, the local Hamiltonian 
is: 

H jk = Ujfij + U k h\ - J jk cos (fa - <f) k ) (22) 

The local approximation to be made here is that, to low- 
est order, the phases on these adjacent sites should be 
locked together. In other words, the degree of freedom 
to be specified is the relative phase <pj — <j) k . This moti- 
vates the introduction of new operators: 

h c = fij + h k 
- _ U k (f>j + Uj<f>k 



Ujhj - U k h k 
Uj + U k 
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j «J k 
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FIG. 3: The link decimation RG step: The crossed link has 
the dominant Josephson coupling. The two sites it joins are 
merged into one cluster, resulting the effective lattice struc- 
ture shown to the right. The cluster C has an effective charg- 
ing energy given in equation (26). 



4>R = 4>j ~ <t>k (23) 
These operators satisfy the commutation relations: 

[<t>c,n c ] = i 



[0. 



(24) 



with all other commutators vanishing. Thus, the trans- 
formation preserves the algebra of number and phase op- 
erators. A subtlety arises for the relative coordinate op- 
erators ur and (J)r because, as defined above, ur need 
not be an integer and <pR G [— 27r, 2tt) as opposed to 
(j)R G [0, 27r). To deal with this difficulty, we may make 
the additional approximation of treating <j)R as a non- 
compact variable. This makes ur continuous instead of 
discrete. Then, in terms of the new cluster and relative 
coordinate operators, the local Hamiltonian (22) reads: 



UjUk 

u* + u k 



h 2 c 



Uk)h 2 R 



Jjk cos ( 



i) (25) 



To lowest order, we set <j)R = 0. This decimation of the 
relative phase leaves the cluster phase <\>c unspecified, so 
two phase degrees of freedom have been reduced to one. 
The first term in Hj k shows that the inverse charging 
energies add like the capacitances of capacitors in parallel 
to give the charging energy for the cluster: 



U c 



1 



UjUk 



J_ 



i 



(26) 



Figure 3 depicts this process of link decimation 13 . 

Higher order corrections to this picture, arising from 
harmonic vibrations of the phases that make up the clus- 
ter, can be obtained by considering the part of the lo- 
cal Hamiltonian involving the relative coordinate. These 
terms act approximately like a simple harmonic oscilla- 
tor Hamiltonian on the basis of ur eigenstates. Thus, 
the ground state for the relative coordinate can be ap- 
proximated by a simple harmonic oscillator ground state: 



\^r) 



7 



4 J-c 



driR exp 



7 s 



\ n R) 



(27) 
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2(Uj + U k ) 



jk 



(28) 



This approximation is used in the numerics to compute 
Debye- Waller factors that modify Josephson couplings 
entering the newly formed cluster. Quantum fluctua- 
tions of (j)R weaken the phase coherence of the cluster, 
and consequently, suppress these Josephson couplings. 
Mathematically, the Debye- Waller factors arise because, 
in writing down the local, two site Hamiltonian (22), we 
have neglected that 4>r also appears in the other links 
penetrating the two sites j and k. Consider a Josephson 
coupling from a third site m to the site j. This corre- 
sponds to a term in the full Hamiltonian (7): 



cos ( 



(j>m ~<t>C - 

(j>m ~ 0c) COS (j1i(/>r) 

in (J> m - C ) sin (jii<f> R 

(j>m - 0c) (COS (jJLi4>j^j 

(j>m ~ 0c)(sin (/ii0 



with: 



(29) 



(30) 



The angle brackets in the final line of equation (29) refer 
to averages taken in the relative coordinate ground state 
(27). The expectation value of the sine vanishes, and the 
expectation value of the cosine yields the Debye- Waller 
factor: 



cdw, 



sin 2 (7r/ij) 



E 



(i 2 



(q 2 



7 s 



(31) 

In the numerics, we truncate the calculation of this sum 
at a specified order, |g m ax| = 20, and multiply the Joseph- 
son coupling J m j by the result to find the new Joseph- 
son coupling J m c penetrating the cluster. Note that the 
Debye- Waller factor for links penetrating the site k is, in 
general, not equal to cdwj, but its calculation is com- 
pletely analogous. 



B. Adaptations for 2D 

Sum Rule 

As shown by Altman, Kafri, Polkovnikov, and Refael, 
the two renormalization steps outlined above form the 
basis of the strong disorder renormalization group for the 
disordered rotor model in ID 13 . In higher dimensions, 
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FIG. 4: Site decimation with the sum rule: The site with 
the dominant charging energy (marked with an X) is coupled 
to two sites (j and k) that are already coupled to one an- 
other. After site decimation, the effective Josephson coupling 
between sites j and k is the sum of the old coupling and the 
effective coupling generated through decimation of site X (see 
equation (32)). 



FIG. 5: Link decimation with the sum rule: The two sites 
connected by the dominant Josephson coupling (sites j and 
k) are both coupled to a third site m. Following link decima- 
tion, the effective Josephson coupling between site m and the 
cluster C is the sum of the two preexisting couplings between 
site m and sites j and k (see equation (33)), up to corrections 
coming rom Debye- Waller factors. 



the geometry of the lattice changes as the renormaliza- 
tion proceeds, and this presents complications. For ex- 
ample, in Figure 4, decimation of a site X produces an 
effective Josephson coupling between two sites j and k 
that are already linked to one another by a preexisting 
Josephson coupling. In our numerics, we sum the pre- 
existing and new coupling to form the effective coupling 
between the remaining sites: 



Jjk — J jk 



JjxJxk 
U X 



(32) 



A similar situation can arise during link decimation. In 
Figure 5, a cluster is formed by two sites j and fc, each of 
which is connected to a third site m. Up to corrections 
coming from Debye- Waller factors, the effective Joseph- 
son coupling joining site m to the cluster is: 



mC 



Jn 



J n 



(33) 



Some authors replace the sum rule in equations (32) 
and (33) with a maximum rule 19 > 26 > 27 . The motivation 
behind the maximum rule is that, in the strong disorder 
limit: 



max 



J 



JjxJxk 



U x 



J 



jk 



JjxJxk 
U x 



(34) 



This should be a good approximation in an infinite dis- 
order context. For our model however, we find that the 
sum rule increases the class of distributions which find 
the unstable fixed point depicted schematically in Figure 
1. For further discussion of the difference between the 
sum and maximum rules, please consult Appendix A. 



Thresholding 

In dimensions greater than one, there is a tendency 
for the numerics to slow down considerably if the renor- 
malization procedure involves a lot of site decimation. 
Again, the source of the problem is the evolution of the 
lattice under the RG. If a site X is decimated, then ef- 
fective links are generated between each pair of sites that 
were previously coupled to X. Thus, site decimation gen- 
erates many new couplings, increasing the coordination 
number of the effective lattice. At the same time, the site 
decimation step takes computer time quadratic in the co- 
ordination number of the site being decimated. To apply 
the procedure to large lattices, it is necessary to find a 
way to circumvent this difficulty. 

At the beginning of the RG, we specify a thresholding 
parameter, which we call a. During a site decimation, if 
a new Josephson coupling is created between sites j and 
k such that: 



Jjk 



JjxJxk 
U x 



< aUx = &Q 



(35) 



then the coupling is thrown away. For convenience in 
implementation, the new bond is ignored only if it does 
not sum with a preexisting Josephson coupling. If a is 
chosen to be very small, then ignoring the coupling will 
hopefully not affect the future course of RG. However, to 
be more careful, it is better to perform an extrapolation 
in the threshold a to see if the numerics converge. Using 
this thresholding procedure, we are able to reach lattices 
up to size 300 x 300, if we additionally require averaging 
over a reasonably large number of disorder samples. In 
this paper, unless otherwise stated, we always use 10 3 
samples for any given choice of distributions. 



Distribution Flows 

Typically, in an application of the strong disorder 
renormalization method, it is interesting to monitor the 
flow of the distributions of the various couplings as the 
RG proceeds. This is straightforward for a ID chain, but 
in higher dimensions, there is yet again a complication 
from the evolving lattice structure. As the renormal- 
ization proceeds, it is possible to generate very highly 
connected lattices. Many of the effective Josephson cou- 
plings will, however, be exceedingly small. Incorporating 
these anomalously small couplings into statistics can be 
misleading. Despite the large number of weak bonds, 
there may exist a number of strong bonds sufficient to 
produce superfluid clusters. In fact, including the weak 
bonds in statistics is analogous to polluting the statistics 
with the inactive next-nearest neighbor Josephson cou- 
plings in the original lattice. It is more appropriate to 
follow Motrunich et al. and focus on the largest O(N) 
Josephson couplings, where N is the number of sites re- 
maining in the effective lattice 27 . In the remainder of the 
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FIG. 6: The projection, in the A J/ J vs. J/U plane, of the 
flows of the coupling distributions at different stages of the 
RG. The initial distributions Pi(U) and Pi (J) are both trun- 
cated Gaussians, and Jo (the center of the initial J distribu- 
tion) is used as the tuning parameter. Each flow corresponds 
to a different choice of the tuning parameter. The flows start 
at the bottom of the figure and go up and to the left or up and 
to the right. A smaller value of the thresholding parameter is 
used near criticality as indicated by the legend. All runs were 
done on L = 100 lattices. 
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FIG. 7: Same as Figure 6, except Pi(U) is Gaussian and 
Pi (J) oc J -1-6 , with cutoffs chosen to make the latter distri- 
bution very broad. The parameter Uo is used to tune through 
the transition. The flows begin near the center of the figure. 
To the left of the figure, flows initially propagate towards the 
bottom left but eventually turn around and propagate to- 
wards the top left. To the right of the figure, flows initially 
propagate towards the bottom left but eventually turn around 
and propagate towards the bottom right. All runs were done 
on L — 300 lattices. 



paper, the "Josephson coupling distribution" will there- 
fore refer solely to the dominant 2N effective Josephson 
couplings at any stage in the RG, and all statistics will 
be done only on these 2N couplings. 



IV. NUMERICAL APPLICATION OF THE 
STRONG DISORDER RG 

In this section, we present numerical data collected 
from our implementation of the strong disorder renor- 
malization group. First, we explore the strong disorder 
RG flows of the distributions of charging energies and 
Josephson couplings. This investigation points to the ex- 
istence of an unstable fixed point of the RG flow. We find 
that the presence of this fixed point is robust to many dif- 
ferent changes in the choices of the initial distributions. 
Next, we examine the distributions generated by the RG 
near this fixed point and find that universal physics arises 
in its vicinity. Subsequently, we proceed away from the 
fixed point to study properties of the phases of the dis- 
ordered rotor model. We find phases that we tentatively 
identify as Mott insulating, glassy, and superfluid, and 
furthermore, we find that the unstable fixed point gov- 
erns the putative glass-superfluid transition. We defer 
detailed interpretation of the data and analysis of the 
transition to Section V. 



A. Flow Diagrams and the Finite Disorder Fixed 
Point 

In Figures 6-8, we plot flows of quantities that charac- 
terize the Josephson coupling and charging energy distri- 
butions. We emphasize again that, in the context of our 
study of distribution flows, the "Josephson coupling dis- 
tribution" actually only includes the greatest 2N Joseph- 




FIG. 8: Same as Figure 6, except Pi(U) oc U 5 and Pi(J) oc 
J -3 . The tuning parameter is J m in, the lower cutoff of Pi(J). 
All runs were done on L — 150 lattices. 



son couplings, where N is the number of sites remain- 
ing in the effective lattice. After M decimation steps of 
the RG, we stop the procedure and look at the remain- 
ing charging energies and these dominant Josephson cou- 
plings. We then use these values to update estimates for 
the mean and standard deviation corresponding to that 
step in the RG. For each realization of the randomness 
(i.e. each sample), we do this for many different choices of 
M, and we repeat the process for 10 3 realizations of the 
randomness. Ultimately, this procedure gives a "flow" 
that characterizes the disorder- averaged evolution of the 
distributions at different stages of the renormalization. 

The x-axes of Figures 6-8 give the ratio of the means 
of the two distributions. Meanwhile, the y-axes give the 
ratio of the standard deviation of the Josephson coupling 
distribution to the mean of the distribution. The plots 
actually show 2D projections of flows that occur in the 
space of all possible distributions. At the very least, these 
plots imply the existence of a third axis, namely 
which may carry important information. Nevertheless, 
these highly simplified 2D pictures are surprisingly ef- 
fective in describing the fate of the model with different 
choices of parameters. In interpreting Figures 6-8, the 
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reader will likely find it helpful to glance back to Figure 
1, which shows a schematization of the flows. 

Figure 6 specifically corresponds to flows for initial dis- 
tributions Pi(U) and Pi(J) that are Gaussian. The cen- 
ter of the Josephson coupling distribution is used as the 
tuning parameter. To the left of the plot are two flows 
that propagate to the top left of the diagram, towards 
small |j and large ~. Since these flows propagate to- 
wards high J7, it is tempting to identify them as flowing 
towards an insulating regime. Meanwhile, to the right 
of the plot, there are seven flows that propagate towards 
high J, and it is tempting to identify these as propagat- 
ing towards a superfluid regime. At the interface between 
these two behaviors, the flows "slow down" and travel a 
shorter distance in the plane. This behavior is suggestive 
of a separatrix flow that terminates at an unstable fixed 
point, as shown in Figure 1. 

Our next goal will be to show that the behavior indi- 
cating the presence of this unstable fixed point is robust 
to changes in the choice of the initial distributions. In 
Figure 7, P%(U) is a Gaussian, and Pi(J) oc J -1,6 . The 
center of the charging energy distribution, Uq, is used 
as the tuning parameter. The numerical choices place 
the flows initially above and to the right of the loca- 
tion of the unstable fixed point in the previous figure. 
From the point of view of the strong disorder RG proce- 
dure, this choice of initial distributions is advantageous, 
because the flows begin in a regime of high disorder in 
J, where the procedure is more accurate. Later in the 
paper, after presenting evidence of the universal physics 
that emerges in the disordered rotor model, we will focus 
on this choice of distributions exclusively. Therefore, we 
have collected additional details about these distributions 
in Appendix C. Note that the leftmost flows in Figure 7 
initially propagate towards the lower left hand corner of 
the figure; then, they turn upward, continuing onward 
to lower |j but now also towards high ~. Hence, they 
share the same qualitative fate as the leftmost flows in 
Figure 6. To the right of Figure 7, the flows initially also 
propagate towards the lower left; however, these flows 
ultimately turn around and propagate towards high ^. 
The separatrix that divides these two classes of flows ap- 
pears to terminate in the same critical region that was 
seen in Figure 6. 

In Figure 8, we make yet another choice of initial dis- 
tributions. Now, Pi(J) oc J' 3 and P^U) oc U 5 . The 
resulting flow diagram again suggests the presence of an 
unstable fixed point in the same critical region. It would 
be misleading, however, to suggest that every flow dia- 
gram generated by the RG will have the nice properties 
of Figures 6-8. We provide a counterexample in Figure 
9, in which Pi(U) is bimodal and P%{J) oc J -16 . Panel 
(a) shows the extremely complicated behavior of some 
of the flows. These features are reflections of the struc- 
tural details of the bimodal distribution. We will see 
shortly that, at least in the vicinity of criticality, the RG 
works to wash away these details and construct univer- 
sal distributions. After these universal distributions are 
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FIG. 9: In these numerical flow diagrams, P%(U) is bimodal 
and Pi (J) oc J -1-6 , with cutoffs chosen to make the latter 
distribution very broad. In panel (a), we show a few sam- 
ple flows that start near the top of the figure and initially 
propagate towards the lower left hand corner. The complex 
features of these flows reflect the structural details of the bi- 
modal U distribution. In panel (b), we exhibit the flows at 
late RG times, when the procedure has had an opportunity 
to renormalize away the microscopic details of the initial dis- 
tributions. Then, the flows are, at least qualitatively, more 
similar to those seen in Figure 7. All runs were done on 
L = 200 lattices. 



somewhat well approximated, the flows should be more 
well behaved, but in Figure 9, we see a non-universal era 
of the flows, where the complexities of the initial distri- 
butions can manifest in complicated flows. To bring out 
this point more clearly, we have removed data for the 
early stages of the RG in panel (b). Now the "late RG- 
time" behavior of the flows falls more nicely in line with 
what is seen in Figure 7. 



B. Universal Distributions 

Near the unstable finite disorder fixed point of the RG 
flow, we expect universal physics to emerge. Certain as- 
pects of the critical behavior should be independent of 
microscopic details, including the structure of the initial 
distributions. The universality of the fixed point should 
become evident in the forms of the renormalized distribu- 
tions generated through the RG: whatever the initial dis- 
tributions may be, they should evolve towards universal 
forms, provided that they put the system near criticality. 

We first focus on determining the universal form of 
the fixed point Josephson coupling distribution. Figure 
10 shows data for the four different choices of the initial 
distributions that we explore in this paper. In panels (a), 
(c), and (d), P%{U) and P%{J) have the same qualitative 
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form, and in panel (b), P^J) ex J -1,6 and P^U) ex U 1 - 6 . 
We tune the parameters characterizing the distributions 
such that the flows propagate near the unstable fixed 
point, run the numerics on 100 x 100 lattices, and plot 
the initial distributions alongside the renormalized dis- 
tributions when 100 sites remain in the effective lattice. 
For the renormalized distributions, we again only include 
the dominant 2N = 200 Josephson couplings for each 
sample. The renormalized distributions suggest that the 
RG indeed washes away the details of the initial choices, 
leaving a power law in each case. The universality of this 
power law is more striking in Figure 11, where we plot 
the renormalized distributions for the four cases together. 
In this plot, we scale J for each of the four cases by the 
mean RG scale Q when only 100 sites remain. This scal- 
ing causes the distributions to nearly collapse onto one 
another, revealing the universal form: 

M£H£)" (36) 

We will momentarily defer providing a numerical esti- 
mate of (/?, in anticipation of presenting higher quality 
data, taken from runs on larger lattices, below. 

We have not plotted the renormalized charging energy 
distributions for the four cases shown in Figure 10. Were 
we to do so, we would see that these distributions, while 
showing hints of universality, are not as strikingly uni- 
versal as the corresponding renormalized Josephson cou- 
pling distributions. The reason for this is the following: 
in three out of the four cases, the initial distributions 
have Jmax < ^max- Several of the initial distributions we 
study in this paper satisfy this property, because in di- 
mensions greater than one, interesting choices of distribu- 
tions typically have most bare charging energies greater 
than most bare Josephson couplings. Otherwise, they al- 
most certainly yield superfluid behavior. Consequently, 
for three out of the four cases in Figure 10, the RG be- 
gins with only site decimations. These site decimations 
dramatically modify the Josephson coupling distribution, 
but the charging energy distribution is, to a large ex- 
tent, only truncated from above by the renormalization 
scale. Later on in the procedure, after many sites have 
been decimated away, the RG enters a regime where the 
charging energy and Josephson coupling scales compete. 
Only then do link decimations begin to occur, and only 
then can the charging energy distribution begin to evolve 
in a nontrivial way. However, by this point, there is far 
less RG time remaining for the fixed point distribution 
to emerge. 

There are two ways to circumvent this difficulty. One 
strategy is to note that this problem of insufficient RG 
time would not arise if we had access to arbitrarily large 
lattices. We could follow the renormalization as long as 
necessary to construct the universal distributions. Thus, 
we can try to explore larger lattices up to the limits set 
by our computational capabilities. On the other hand, 
another solution is to work with very wide initial distri- 
butions of Josephson couplings. These are distributions 




ln(J) 

FIG. 10: Log- log plots of initial and renormalized Josephson 
coupling distributions for near-critical flows. All runs were 
done onL= 100 lattices with a = 5x 10 -6 . Each plot shows 
the initial distribution and the distribution when the effec- 
tive lattice has 1% of the original number of sites. The initial 
distributions have four different forms, but the distributions 
after renormalization show a universal power law. Note that 
the plots of initial distributions in these plots were not con- 
structed from actual data (i.e. actual numerical sampling of 
the distributions) but were instead constructed by hand. 
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FIG. 11: The distributions from Figure 10, with the Joseph- 
son coupling strength scaled by the mean RG scale Q at the 
corresponding stage of the RG. The near collapse of the distri- 
butions reveals the universal power law form of the Josephson 
coupling distribution near the unstable fixed point. A refined 
version of this plot, showing data for larger lattices (but also 
different choices of initial distributions) , can be found in panel 
(b) of Figure 13. 



which have large ~. As such, they correspond to flows 
that begin above the unstable fixed point in our ~ vs. 

jj flow diagrams. Using such distributions, it is possible 
to engineer situations where most bare charging energies 
exceed most bare Josephson couplings but where, due 
to the presence of a small fraction of anomalously large 
Josephson couplings, J max > ^max at the beginning of 
the RG. If the parameters are chosen appropriately, the 
renormalization procedure will begin with a few link dec- 
imations, until the charging energy and Josephson cou- 
pling scales meet. After this point, the RG will feature 
an interplay of site and link decimations. Thus, both 
the Josephson coupling and charging energy distributions 
will evolve nontrivially. 

To target the fixed point charging energy distribution, 
we apply both of the strategies. We proceed to 200 x 200 
lattices and compute renormalized distributions when the 
effective lattice has 200 sites remaining. Additionally, we 
work with very wide initial Josephson coupling distribu- 
tions. In order to achieve large initial ~, we restrict our 
attention to power law distributions of Josephson cou- 
plings of the form Pi(J) oc J -1-6 . We vary the choice of 
the initial charging energy distribution and tune the pa- 
rameters near criticality. The results are shown in Figure 
12. Now, the RG does generate a universal form for the 
charging energy distribution, and in Figure 13, we scale 
the renormalized distributions by the corresponding RG 
scales to expose the universality more clearly. Figure 
13 suggests that the functional form of the fixed point 
charging energy distribution may be: 



oc 



exp 



m >< (§) 



(37) 



where gu is an energy scale below which the charging 
energies are exponentially rare. We have been unable to 
extract a good estimate of /3. Panel (a) of Figure 13 pre- 
supposes ft « 1, and the linearity of the plots suggests 



that this may be close to the actual value. Taking (3 = 1 
and focusing on the case where Pi(U) is Gaussian (be- 
cause that is the choice of initial distributions for which 
we have most accurately targeted criticality), we fit: 



9u_ 



0.66 ±0.02 



(38) 



We should note that qualitatively similar charging energy 
distributions to those seen in Figure 12 still emerge near 
criticality if we relax the restrictions of initially power law 
J distributions and initially high ~ . This is true of the 
distributions studied in Figure 10, even in the fiat and 
bimodal cases where J ma x < C/min initially and clusters 
can only form due to the use of the sum rule. As argued 
above, the additional restrictions we impose on Pi(J) in 
Figure 12 simply allow the RG to construct the fixed 
point distributions more cleanly. 

In the lower panel of Figure 13, we additionally present 
data for the Josephson coupling distributions at the same 
stage of the RG. Again focusing on the data for the case 
where P%(U) is Gaussian, we can estimate: 



(p : 



1.15 ±0.01 



(39) 



Before proceeding, we should note that the form of the 
fixed point J distribution allows us to construct an ar- 
gument for the validity of the RG near criticality. We 
expand upon this argument greatly in Appendix D, but 
we will sketch the basic premise here. Essentially, we 
should consider the implications of the fixed point distri- 
butions for the reliability of each of the RG steps. The 
validity of site decimation rests on the reliability of the 
perturbative treatment of the Josephson couplings en- 
tering the site with the dominant charging energy. The 
form of the critical Josephson coupling distribution im- 
mediately guarantees that most Josephson couplings are 
much weaker than the RG scale. For the Gaussian case 
in Figure 13, the ratio of the median J to the RG scale 
is ^0.11±0.01. Hence, the site decimation is usu- 
ally very safe. In the case of link decimation, a simi- 
lar argument allows us to ignore, to leading order, other 
Josephson couplings penetrating the sites joined by the 
dominant coupling. However, the structure of the crit- 
ical charging energy distribution, shown in Figure 13, 
actually suggests that there can be a large number of 
charging energies of the same order as the RG scale; in 
particular, the ratio of the median U to the RG scale 
is « 0.67 ± 0.01. Consequently, the question of 

the reliability of link decimation reduces to the follow- 
ing: in a single junction (or two-site) problem, how re- 
liable is cluster formation when both charging energies 
are weaker than the Josephson coupling but potentially 
of the same order-of-magnitude? We address this ques- 
tion in Appendix D and find that the link decimation 
step also seems to be reasonably safe. 
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FIG. 12: Log- log plots of initial and renormalized charging en- 
ergy distributions for near-critical flows. All runs were done 
on L — 200 lattices with a = 5 x 10 -6 . Each plot shows 
the initial distribution and the distribution when the effective 
lattice has 0.5% of the original number of sites. The initial 
charging energy distributions have four different forms, but 
the distributions after renormalization show a universal form. 
Note that the plots of initial distributions in these plots were 
not constructed from actual data (i.e. actual numerical sam- 
pling of the distributions) but were instead constructed by 
hand. 
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FIG. 13: In panel (a), the renormalized distributions from 
Figure 12 are plotted together, with the charging energies U 
scaled by the mean RG scale Q. In panel (b), we provide a 
similar plot for the renormalized Josephson coupling distribu- 
tions produced by the runs in Figure 12. 



C. Physical Properties and Finite Size Scaling 

To determine a preliminary classification of the phases 
of the model, we now measure four physical properties. 
First, we measure s max , the size of the largest cluster 
formed by link decimations during the renormalization 
procedure. This corresponds to the largest domain of 
superfluid ordering in the system. We also measure 52, 
the size of the second largest cluster. We will see that the 
behavior of this quantity differs dramatically from that 
of 

5 max i n the superfluid phase, and therefore, both are 
interesting quantities to measure. 

We also calculate the charging gap for the system, 
A m i n . We can estimate this quantity by simply mea- 
suring the charging energy of the final site remaining in 
theRG. 

Finally, we measure a susceptibility towards superfluid 
ordering. Consider adding a perturbation of the following 
form to the rotor model Hamiltonian: 



H' = -hJ2 cos (h) 



(40) 



In the RG, we can evaluate the linear response of the sys- 
tem to this perturbation and measure the susceptibility: 



X 



1 ^ d (cos (j)j) 

J2 1^ _ 



dh 



(41) 



/i=0 



The terms of this sum are computed during site decima- 
tion. Perturbation theory gives a single site susceptibility 
of jj^, where X is the site being decimated. Neglecting 
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corrections from harmonic fluctuations, we can find the 
contribution of a cluster to the susceptibility by multi- 
plying the perturbative result by s 2 , where s is the size 
of the cluster. One factor of s arises from the fact that 
the cluster represents s terms in the original sum (41), 
and the other follows from the fact that the effective field 
coupling to the cluster phase is enhanced s times when s 
phases rotate together. When harmonic fluctuations are 
taken into account, both of these factors of s should be 
replaced by a renormalized factor which we denote as b. 
This 6-factor accounts for the fact that quantum fluctu- 
ations weaken the phase coherence of the cluster. For a 
bare site, 6 = 1, and when two sites j and k merge, the 
renormalized 6- factor for the cluster C is: 



(a) 



be = bjCnwj + bkCD\v,k 



(42) 



where cdwj an d cnw,k are the Debye- Waller factors (31). 
Hence, the total contribution of the cluster to the suscep- 
tibility, before the normalization by , is: 



Xc 



U c 



(43) 



where be and Uc are the 6- factor and charging energy of 
the cluster respectively. Further details of this calcula- 
tion can be found in Appendix B. 

For each lattice size, we obtain an estimate for the four 
quantities s max , 52, A min , and \- Then, we examine how 
these estimates vary as we raise L. For certain types of 
distributions, computational limitations force us to work 
on relatively small lattices. This occurs, for example, 
when both Pi(U) and Pi(J) are bimodal, and we study 
this case in Figure 14. 

In panel (a) of Figure 14, there is no cluster formation 
whatsoever. Hence, s max = S2 = 1. This results in a gap 
A m i n that remains asymptotically constant as it cannot 
be lower than the lower bound of the initial charging 
energy distribution. The susceptibility \ also remains 
asymptotically constant. 

Next, in panel (b), we find a regime in which link dec- 
imations do occur and clusters do form. Moreover, s max 
and S2 grow with system size, with what appears to be 
a power law for the largest lattice sizes that we explore. 
This growth is, however, sub quadratic in L, meaning that 
s max grows subextensively with lattice size. Meanwhile, 
A m i n decays with a power slower than L -2 , and the sus- 
ceptibility x remains constant with L. 

In panel (c), all quantities show power law behavior out 
to L = 100, including the susceptibility which appears to 
grow with a very slow power. The growth of s max is still 
slower than L 2 , so the largest cluster is still subextensive. 

Finally, in panel (d), we find a regime in which s max 
grows as L 2 , reflecting the formation of macroscopic clus- 
ters that scale extensively with the size of the lattice. The 
gap A min decays as L -2 , and the susceptibility shows an 
approximately L 4 growth for large L. Perhaps surpris- 
ingly, although 82 grows with system size, it does so more 
slowly than it does in panel (c). 
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FIG. 14: Four characteristic behaviors of s ma x, $2, A m i n , and 
X as a function of system size. Here, P%(U) and P%(J) are 
both bimodal, and Jh, the center of the higher peak of the 
Josephson coupling distribution, is used as the tuning param- 
eter. All quantities have been normalized by their value for 
L — 10, and data is shown for L — 10 to L — 100. In panel 
(a), the quantities reflect the purely local physics of the Mott 
insulator. In panel (b), s max and S2 grow subextensively with 
system size, with what appears to be a power law. The gap 
also decays with a slow power, and the susceptibility remains 
constant. In panel (c), all quantities show power law behavior 
out to L = 100. The reference line shows L 2 growth. Panel 
(d) reflects the macroscopic clusters of the superfluid phase. 
The cluster size s max is parallel to the L 2 reference line, and 
the susceptibility % is parallel to the L 4 reference line for large 
L. 



We now turn to a class of distributions for which we 
can reach larger lattice sizes. In particular, we return to 
the data set described in Appendix C, in which Pi(U) 
is Gaussian and Pi(J) oc J~ 16 . Data for this choice of 
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distributions is shown in Figure 15. 

Panel (a) of this figure qualitatively reproduces the 
features of panel (a) of Figure 14. Panel (b) of the new 
figure, on the other hand, differs from panel (b) of Figure 
14 in an important way. For large L, the power law 
behaviors of s max , S2 and A m i n are lost, and all three 
quantities vary more slowly. In Figure 14, this effect may 
have been hidden by the use of smaller system sizes. 

If the parameters are tuned such that the correspond- 
ing flow propagates very close to the unstable fixed point, 
then we do find a regime in which all quantities show 
nearly power law behavior out to L = 300. This regime 
is depicted in panel (c) of Figure 15. 

Tuning past this point, we enter a regime in which 
macroscopic clusters form. Panel (d) of Figure 15 shows 
the behavior of the four quantities in this regime, and we 
see that most of the essential features of the correspond- 
ing panel of Figure 14 are reproduced. An important 
feature of the plot in panel (d) is that we can clearly see 
that the behavior of S2 in this regime is closer to that 
observed in panel (b) than in the intervening panel (c). 

The plots in Figures 14 and 15 are suggestively labeled 
with their corresponding phase identifications. We will 
provisionally use these labels for convenience in referring 
to these regimes, in advance of presenting arguments for 
these classifications in Section V. In the flow diagrams 
that we presented earlier, choices of initial distributions 
that correspond to the Mott insulating and glassy behav- 
iors shown in panels (a) and (b) of Figures 14 and 15 flow 
to the stable insulating region where jj — » 0. The super- 
fluid behavior in panel (d) of the figures corresponds to 
flows towards the high ^ regime of the flow diagrams. 
The pure power law behavior of panel (c) emerges for 
flows that propagate very close to the proposed unstable 
fixed point. This suggests that this fixed point may con- 
trol the glass to superfluid transition of the disordered 
rotor model. 

For now, we assume that this is the case and investi- 
gate more closely the behavior of physical quantities in 
the vicinity of this proposed transition. Having provided 
evidence of the universality that emerges near the critical 
point, here and in the remainder of this paper, we will 
focus exclusively on the choice of distributions detailed 
in Appendix C. In Figure 16, we show that plots of phys- 
ical quantities vs. tuning parameter, taken for different 
L, can be collapsed onto universal curves. We will use 
this scaling collapse in Section V to determine the criti- 
cal exponents governing the putative transition between 
glassy and superfluid phases. 



D. Cluster Densities and 6-factors 
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FIG. 15: Four characteristic behaviors of s ma x, S2, A m i n , and 
X as a function of system size. The initial distributions are 
those described in Appendix C. All quantities have been nor- 
malized by their value for L — 25, and data is shown for 
L = 25 to L = 300. In the four panels, U = 400, 9.2, 8.97, 
and 8.8 respectively. Panel (a) reflects the purely local physics 
of the Mott insulator. Panel (b) shows a glassy regime char- 
acterized by rare-regions clusters that grow subextensively in 
system size. The reference line shows the power law that s max 
obeys near criticality. This nearly critical regime is shown in 
in panel (c). The reference line here shows L 2 growth. Fi- 
nally, panel (d) shows the superfluid regime, in which s max is 
asymptotically parallel to the L 2 reference line. The suscep- 
tibility x is expected to grow as L 4 for very large L, but it 
does not quite reach this behavior (indicated with a reference 
line) for L < 300. 



The physical property that distinguishes the Mott glass 
from the Bose glass is the compressibility. Later, we will 
show that, in order to calculate the compressibility, it is 
insufficient to consider just the minimum charging gap. 
Within each sample, the RG may form several clusters, 



each of which implies a local gap for adding and removing 
bosons. We will need to monitor all of these local gaps 
to find the compressibility. 

More specifically, in this section, we will calculate the 
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FIG. 17: The number per unit area of clusters with local gap 
near A for a choice of parameters in the glassy phase. The 
initial distributions are those described in Appendix C, with 
Pi(U) Gaussian and P%(J) oc J -1-6 . The value of the tuning 
parameter is Uo — 8.8, and this puts the system on the glassy 
side of the transition at Uo ~ 8.97. Data is shown for L — 75, 
150, and 300 lattices. The density decays faster than a power 
law at small A. 



FIG. 16: Panel (a) shows the scaling collapse of Smax as a 
function of tuning parameter, and panel (b) shows a similar 
collapse of S2- Each line corresponds to a different value of 
the lattice size. We show data for L = 25, 50, 75, 100, 150, 
200, and 300. The insets show magnified views of the vicinity 
of the critical point for the four largest lattice sizes. The error 
bars, which are difficult to see in the larger plots, are clearly 
visible in the insets. To see the values of the exponents v and 
df used in each panel, consult equations (71), (72), (76), and 
(77). 



density (number per unit area) of clusters of a given size 
and of clusters with gaps in a given range of energy. We 
call these densities p(s) and p(A) respectively. The latter 
quantity gives a "density-of-states" for addition of single 
bosons or holes to the system. For a single sample cor- 
responding to a specific choice of initial distributions, we 
monitor the size and local charging gap for all the clus- 
ters formed during the renormalization, excluding bare 
sites that are not clustered by the RG. In the case of the 
local charging gap, we again estimate this quantity as the 
charging energy of a cluster at the time of decimation. In 
principle, we could also include perturbative corrections 
to this local gap, but we omit these and do not anticipate 
that they would affect the behavior of the density at low 
A. We pool data for 10 3 samples, choose a discretization 
of energy to determine a histogram bin size, compute a 
histogram of gaps and a histogram of cluster sizes, and 
finally normalize these histograms by the total simulated 
surface area: L 2 times the number of samples. 

Our study of these densities will bring into focus the 
crucial difference between two types of clusters formed 
by the RG: rare-regions clusters and the macroscopic 
clusters that characterize the superfluid phase. We will, 
therefore, also take the opportunity to examine how the 
b- factors, which quantify the effect of harmonic fluctua- 
tions on the susceptibility, vary as a function of s for the 
two types of clusters. 



The Charging Gap Density p(A) 

Note that p(s) and p(A) are not particularly inter- 
esting for choices of distributions and parameters that 
yield the Mott insulating behavior from Figures 14 and 
15. The profile of p(A) will be identical to the profile 
of the initial charging energy distribution, and because 
we choose this distribution to be bounded from below by 
some positive J7 m i n , it can be shown that this always cor- 
responds to an incompressible phase. Hence, we begin 
by focusing on the glassy regime. 

Figure 17 shows the gap density for a choice of pa- 
rameters in the glassy phase. As the size of the lattice 
is raised, the density profile remains essentially invari- 
ant at large A, but smaller gaps, corresponding to larger 
clusters, begin to appear. However, these smaller gaps 
simply fill out a decay to as A — >• 0. 

Now, we turn our attention to the putative superfluid 
phase. The gap density in this phase is shown in Figure 
18. In panel (a), there is an invariant piece to the dis- 
tribution, but at very low A, a second peak appears as 
well. Panel (b) of Figure 18 shows a magnified view of 
this low A peak. This peak propagates towards A = as 
the system size is raised. Accompanying the propagation 
is a shrinking: the integrated weight of the low A peak 
shrinks as L~ 2 . The consequences of these two effects 
need to be taken into account carefully to calculate the 
compressibility. 

The Cluster Size Density p(s) 

We now consider how p(s), the density of clusters of 
size 5, varies as we sweep through the glassy regime and 
into the superfluid. Panel (a) of Figure 19 shows the ap- 
proach to the transition from the glassy side. Very close 
to the transition at Uo ~ 8.97, p(s) exhibits a striking 
power law decay. Proceeding into the proposed glassy 
phase, the power law decay persists at small s. However, 
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FIG. 18: In panel (a), the density (number per unit area) of 
clusters with a given gap A for a choice of parameters in the 
superfluid phase. The initial distributions are those described 
in Appendix C, with P%{U) Gaussian and P%(J) oc J -1-6 . The 
value of the tuning parameter is Uo — 9.2, and this puts the 
system on the superfluid side of the transition at Uo ~ 8.97. 
Data is shown for L = 75, 150, and 300 lattices. The density 
profiles exhibit two peaks. The broad peak that is visible in 
panel (a) remains invariant with L. To expose the second 
peak, we provide a magnified view of the low A part of the 
density in panel (b). This peak simultaneously shrinks and 
propagates towards A = as the system size is raised. 



FIG. 19: Sweeps of p(s) as the system is tuned through the 
superfluid-insulator transition on L — 300 lattices. The initial 
distributions are those described in Appendix C, with P%(JJ) 
Gaussian and Pi(J) oc J -1 ' 6 . The tuning parameter is [To, 
the center of Pi(U). Panel (a) shows the sweep from deep 
in the glassy phase (Uo = 20) to the transition (Uo = 8.97). 
The closest data set to the transition is indicated by the black 
line with large data point markers. This line is reproduced in 
panel (b) , which shows the sweep from the transition into the 
superfluid phase (up to Uo = 6.5). In the superfluid phase, 
the density plot is characterized by a peak at large s and a 
remnant decay at low s. 



this behavior is cut off by some scale s, beyond which 
p(s) decays very rapidly, essentially exponentially. 

In panel (b), we proceed in the other direction from 
the putative transition, into the superfluid regime. A 
peak, corresponding to the macroscopic clusters, appears 
at large s. The macroscopic cluster in each sample is 
dressed by rare-regions clusters, and these clusters are 
represented by the remnant decay at small s. While the 
large s peak is related to the low A peak in Figure 18, 
the remnant decay at small s is the analogue of the high 
A feature that stays invariant with system size. The low 
s decay in panel (b) of Figure 19 qualitatively resembles 
the decay well inside the glassy regime. In summary, p(s) 
exhibits a pure power law decay in the vicinity of the pro- 
posed glass-superfluid transition; tuning away from crit- 
icality in either direction, and excluding the macroscopic 
clusters of the superfluid phase, the power law form of 
p(s) only survives up to a scale s. For s > s, clusters 
become exponentially rare. 

A type of scaling collapse can be performed for the 
p(s) curves from Figure 19, and this collapse is shown in 
Figure 20. We will see below that this collapse gives a 
complementary set of critical exponents which are related 
by scaling to those that can be extracted from Figure 16. 








30 40 



FIG. 20: Scaling collapse of the p(s) curves from Figure 19. 
Small cluster sizes (s < 30) need to be discarded, because 
they are non-universal. Large cluster sizes (s > 100) need 
to be discarded, because they are noisy. Then, the remaining 
p(s) curves, taken for different values of the tuning parameter, 
collapse onto a universal curve. 



17 



Susceptibility b-factors 

The data presented above for the cluster densities 
p(s) and p(A) highlights the difference between the rare- 
regions and macroscopic clusters generated by the RG. A 
study of how the 6-factors for the clusters vary as a func- 
tion of s can bring out another difference between the 
two types of clusters. Recall that the b- factor was intro- 
duced in equation (43) to quantity the effect of harmonic 
fluctuations on the susceptibility of a superfluid cluster. 
As such, understanding the behavior of the 6-factors will 
also be essential in explaining the behavior of x m the 
various phases of our model. 

In Figure 21, we plot the average value of b for a cluster 
of size s and plot it against s. Again, we work with 
L = 300 lattices and the choice of distributions described 
in Appendix C. Panel (a) shows data for the glassy phase 
and for the non-macroscopic clusters of the superfluid 
phase. We see that b varies with s as a power law: 



(a) 



b(s) ~ 



(44) 



and that the power is consistent for seven different choices 
of the tuning parameter /7 . We will provide an estimate 
of ( in Section V. Panel (b) of Figure 21 shows data for 
the macroscopic clusters when Uq = 8.8. Now, b(s) oc s. 
This behavior can be anticipated from a simple picture 
of macroscopic cluster growth, which we will discuss in 
Section V. 



V. THE PHASES AND QUANTUM PHASE 
TRANSITIONS OF THE DISORDERED ROTOR 
MODEL 

Having collected representative numerical data in the 
previous section, we now assess what the data tells us 
about the phases and phase transitions of our model. Our 
first task will be to confirm the association of phases with 
the behaviors of physical properties that we observed in 
Figures 14 and 15. To this end, we will have to pre- 
emptively introduce one of the main conclusions of the 
discussion below: that our data points to a superfluid- 
insulator transition driven by a percolation- type process. 
The rare-regions clusters of the glassy phase find one an- 
other, and their phases cohere, producing a macroscopic 
cluster of superfluid ordering and driving the transition 
to long range order and global superfluidity. 

Motivated by this intuitive picture of the transition, 
the logic of the discussion below will be the following: 
in the proposed glassy and superfluid regimes of Figure 
15, the cluster size density p(s) exhibits the universal 
features that we would expect from non-percolating and 
percolating phases of standard models of percolation. We 
can use these correspondences to extrapolate the behav- 
iors seen in Figure 15 to the thermodynamic limit, in 
the process showing that these phases will indeed have 
the properties expected of glassy and superfluid phases. 
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FIG. 21: The behavior of the mean b- factor for clusters of size 
s as a function of s. In panel (a), we show data for the glassy 
regime and for the non-macroscopic clusters of the superfluid 
regime. The initial distributions are those described in Ap- 
pendix C, and data is shown for seven values of the tuning 
parameter Uo on L = 300 lattices. The log- log plot shows 
power law behavior of b(s) vs. s. In panel (b), we show data 
for the macroscopic clusters when Uo = 8.8. The plot indi- 
cates that b(s) rsj s for macroscopic clusters. In both plots, the 
error bars indicate the error on the mean b(s) over all clusters 
of size s. For some of the largest and smallest values of s in 
each plot, the absence of an error bar indicates that only one 
cluster of that size was generated in all of the samples. 



Furthermore, we can analyze the critical point and ex- 
tract the critical exponents that govern the superfluid- 
insulator transition. After characterizing this transition, 
we will finally return to the question of the identity of 
the glassy phase and determine if a Mott glass is present 
in the disordered rotor model (7). 

Before proceeding, we should clarify that, although our 
RG produces clusters with size distributions similar to 
models of percolation, our transition is not the standard 
percolation transition. Indeed, the exponents that we re- 
cover are significantly different from the percolation ex- 
ponents on a 2D square lattice. However, the analogy 
to percolation allows us to easily identify the relations 
between the various exponents and the scaling functions 
for the observables. 



A. Phases of the Model 

Mott Insulator 

We briefly digress to describe the phase of our model 
in which the percolation picture is irrelevant, simply be- 
cause there are no regions of superfluid ordering. In a 
Mott insulating phase, the system wants to pin the num- 
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ber fluctuation to zero on each site to avoid the energetic 
costs of charging. Hence, s max and S2 simply stay pinned 
at one as L is raised. Meanwhile, A m i n should asymptote 
to a constant, reflecting the fact that the Mott insulator 
is gapped. A phase without any cluster formation can be 
described by completely local physics. This means that 
the susceptibility can be approximated by disorder aver- 
aging a single site problem. In other words, x should also 
stay constant as the system size is increased. Thus, in 
a Mott insulating phase, we expect exactly the behavior 
seen in panel (a) of Figures 14 and 15. 



Glass 

In non-percolating phases of standard models of per- 
colation, the cluster size density is expected to go as: 



P(s) 



cs 



(45) 



where c is a constant. The function f(x) is expected to 
be approximately constant for x < 1 and to rapidly decay 
for x > l 28 . Equation (45) is consistent with what we 
have observed in our proposed glassy phase in panel (a) of 
Figure 19 and is also consistent with the expectation that, 
in a Griffiths phase, the frequency of occurrence of large 
rare-regions decays exponentially in their size 25 . In our 
numerics, f(x) seems to follow a pure exponential decay 
f(x) ~ e~ x , but the conclusions we present below would 
be qualitatively unchanged if, for example, f(x) ~ e~ x . 

With the form of equation (45) in hand, we can now 
formulate a simple argument for the asymptotic behavior 
of 5 max . In particular, the order-of-magnitude of s max is 
set by the condition 28 : 



L 

E 



Pl(s) « 1 



(46) 



If the left hand side of equation (46) is less than 1, then 
it is unlikely that even one cluster of size greater than 
or equal to s max will be present in a sample of size L 2 . 
In equation (46), Pl($) is the finite size approximant to 
the function p(s) that appears in equation (45). The up- 
per limit of the sum in equation (46) is taken as L 2 be- 
cause larger clusters simply cannot appear in the finite 
size sample. With enough sampling of the random dis- 
tributions, it should, in principle, be able to accurately 
represent the approximant pl(s) out to s = L 2 . The 
data indicates that Pl(s) will simply reproduce p(s) out 
to this value of s, so in the calculations below, we can 
replace the approximant pl(s) by p(s). This will not be 
possible in the superfluid phase. 

For large L, where we also expect s max > s, we can 
use equation (46) to find s max by computing 



L 2 £ p(s) « V 
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J Syr 



dsp(s) 



L 

cL 2 / dsexp 
= csL 2 



cL z I dsexp 

L : 
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r ln(s) 



s 



(47) 



Setting this expression equal to 1 and inverting for «s max , 
we find that, asymptotically in L: 



InL 



(48) 



For large clusters, the link decimation rule for addition of 
charging energies (26) implies that U ~ <s _1 , and there- 
fore: 



1 

InL 



(49) 



An entirely analogous condition to equation (46) can 
be formulated for «S2- We simply replace the right hand 
side of equation (46) with two, indicating that we want 
to find the value of s such that there are likely to be 
two clusters of that size or greater in a typical sample. 
However, the remainder of the calculation is qualitatively 
unaffected by this change. Therefore, 52 should also grow 
logarithmically in this regime: 



S2 



InL 



(50) 



Finally, we turn to the susceptibility. This can be com- 
puted as follows: 



c 

- 2 J2 p(s)L 2 (b(s)) 2 s 

S=l 

•Smax 

X>(s)(K*)) 2 s 

6 = 1 



(51) 



In this calculation, the sum over C is a sum over clus- 
ters, with sc being the size of the cluster. Then, the sum 
over s is, as before, a sum over cluster sizes, and b(s) is 
the average value of the 6-factor for a cluster of size s. 
Figure 21 shows that, in the glassy regime, b varies as a 
power of s and that this power remains the same all the 
way up to criticality. While we do not have a complete 
understanding of this behavior, we can still understand 
the asymptotic behavior of x by reasoning that b can, at 
most, grow linearly in s. This follows from an interpre- 
tation of the 6-factor as the effective number of rotors 
that order with the field in a cluster of size s. Since p(s) 
decays exponentially for large s while b(s) grows at most 
as a power, the sum (51) converges, and x should be 
asymptotically constant: 



lim x = Xo 

L— )-oo 



(52) 
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All of these behaviors are consistent with what has 
been observed numerically in panel (b) of Figure 15. 
Moreover, since logarithmic behavior can be difficult to 
discriminate from a slow power law at low L, they are 
also consistent with panel (b) of Figure 14. Thus, under 
the numerically justified assumption that this regime cor- 
responds to a non-percolating phase, we can deduce that, 
as L increases, arbitrarily large rare-regions of superfluid 
ordering will appear, driving the gap to zero. However, 
the typical size of these regions grows extremely slowly 
(i.e. logarithmically) with system size. The fraction of 
sites occupied by the largest cluster in a typical sample 
vanishes as L — >> oo, so there is no long range order. The 
behavior of this phase for large L corresponds to what 
we should expect for a glassy phase. 



Critical Region 

At the critical point of a percolation transition, the 
cluster size scale in equation (45) is expected to diverge 
as: 



s~\g-g c 



(53) 



where g is the tuning parameter for the transition. This 
divergence is related to the divergence of a correlation 
length which indicates the typical linear extent of the 
largest clusters: 

Z~\g-9cr (54) 

By equation (53), p(s) is a power law at criticality: 

p(s) = cs~ T (55) 

This means that the critical point is characterized by a 
scale invariant, fractal structure of clusters 28 . In turn, 
this implies that £ and s are related by a fractal dimen- 
sion: 

s = f d / (56) 
Equations (53), (54), and (56) together imply: 

1 



(57) 



We will use this scaling law in our analysis of the transi- 
tion below 28 . 

For the present purposes, note that equation (55) is 
once again consistent with what we have observed nu- 
merically in Figure 19. Now, the calculation for s max 
becomes: 

L 2 P( s ) ~ cl2 / dss ~ T = 1 ( 58 ) 

which, when inverted, yields: 

In s max = In L In 



(59) 



r — 1 r — 1 \ c 

c 



In 1 + 



(r - 1)L 2 (-- 2 ) 



Asymptotically, as long as r > 2, this simply corresponds 
to a power law growth of <s max : 



(60) 



On the other hand, since df is the exponent that connects 
length and cluster size scales at the transition, equation 
(60) yields another scaling relation: 



df 



1 



(61) 



Equations (57) and (61) are the usual scaling laws con- 
necting exponents at a percolation transiiton 28 . 

In accordance with equation 60, the gap should close 
as: 



A min ~ L T-i 



= L~ df 



(62) 



Furthermore, as in the glassy regime, the qualitative be- 
havior of the second largest cluster size 52 should be iden- 
tical to that of 5 max : 



$2 ~ L T ~ 1 



(63) 



Turning to the susceptibility, we find that it no longer 
converges to a finite value. At criticality, power law 
behavior of b(s) follows naturally from scale invariance. 
When p(s) ~ s~ T and b(s) ~ s^: 



X 



s=l 



~ s 



dss 1 ^' 

2+2C-T 



max 

4 + 4<C-2r 



(64) 



~ L 

= L d f (1+20-2 

From Figure 21, we can estimate the exponent £: 

0.68 ±0.01 (65) 

With another choice of initial distributions (P%(U) bi- 
modal and P^J) ex J" 16 ), we have found a similar value 
for £. If df(l + 2Q > 2, then x asymptotically diverges 
as a power law, as seen in Figure 15. We will provide 
an estimate of df shortly in equation (72). For now, we 
note that the power observed for x vs - L m Figure 15 is 
consistent with this estimate of df and the estimate for 
( that is given above. In Figure 14, the relatively small 
system sizes likely put us out of the scaling regime for x, 
and this is probably responsible for the extremely slow 
growth of x w ^h L. 



Superfluid 

The percolating phase is characterized by the presence 
of a macroscopic cluster that grows with the size of the 
lattice, so trivially: 



cx L 



(66) 
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and therefore: 



oc L 



(67) 



It is possible to construct an argument along the lines of 
condition (46) for the behavior of s max in equation (66), 
but in this case, it is important not to substitute the 
infinite lattice density p(s) for the finite size approximant 
Pl(s). The subtlety lies in the high s peak observed in the 
density plots in panel (b) of Figure 19. Consistent with 
their low A counterparts in the plots of the gap density, 
the location of these peaks propagates as L 2 towards high 
s as L is raised. Simultaneously, the integrated weight of 
the peaks shrinks as L -2 , reflecting the fact that there is 
only one macroscopic cluster in each sample. Balancing 
the shrinking and propagation, it is possible to see that, 
in order to achieve the condition (46), s max must scale as 
shown above in equation (66). 

The reasoning above has important implications for 
the behavior of 52- Since the weight of the high s peak 
of Pl( s ) shrinks as L -2 , the second largest cluster must 
be drawn from the remnant low s decay. The physical 
picture behind this low s decay is the following: suppose 
we remove the sites belonging to the macroscopic cluster 
from the original lattice. In doing so, we remove some of 
the lowest charging energies and highest Josephson cou- 
plings from the bare lattice. Consequently, the remnant 
lattice is globally insulating. Nevertheless, rare-regions 
of superfluid ordering can form exactly as in the glassy 
phase. It follows that S2 will grow with L as in the glassy 
phase: 



S2 ~ In L 



(68) 



This is responsible for the peak in S2 at criticality. 

To calculate the susceptibility, we first take into ac- 
count the contribution of the macroscopic cluster. The 
behavior of b(s) for a macroscopic cluster can be in- 
ferred from a simplified picture of a large cluster merging 
with single neighbors. As the macroscopic cluster grows, 
its charging energy becomes smaller, driving the Debye- 
Waller factor for the cluster to one. At the same time, the 
cluster's Josephson couplings to other sites grow through 
link addition processes. This means that the Debye- 
Waller factor for another site that is merging with the 
macroscopic cluster also approaches unity. Therefore, 
the 6-factor addition rule (42) approximately becomes 
be = be + 1, and be ~ sc follows. Then: 



Xn 



1 x (&(w)) 2 

L 2 Amir, 



L 2 X A n 



(69) 



The rare-regions clusters dressing the macroscopic cluster 
add a subleading contribution to the susceptibility, which 
we call Xrr- The reasoning of equation (52) indicates 
that this contribution should be asymptotically constant. 





Mott Insulator 


Glass 


Critical 


Superfluid 


Smax 


const. 


ln(L) 


L d f 


L 2 


S2 


const. 


ln(L) 


L d f 


ln(L) 


A m i n 


const. 


1 

ln(L) 


L~ d f 


L- 2 


X 


const. 


const. 


^(1+20-2 


L 4 



TABLE I: Large lattice size (L) behavior of physical proper- 
ties in the three phases and at the critical point of the glass- 
superfluid transition. 



Thus, the quartic behavior of equation (69) is the correct 
asymptotic behavior. Finite size corrections from \rr will 
modify this behavior, however, and this is probably why 
we do not quite see x reach the L 4 behavior in panel (d) 
of Figure 15. 

We have now provided arguments for the behaviors 
observed in each panel of that figure, and we summarize 
this information in Table I. 



B. Quantum Phase Transitons 

Glass- Superfluid Transition 

Earlier, we remarked in passing that systems exhibit- 
ing the behaviors that we have now identified as Mott 
insulating and glassy eventually propagate towards the 
putative insulating region to the top left of the numeri- 
cal flow diagrams. Correspondingly, the systems exhibit- 
ing superfluid behavior propagate towards the putative 
superfluid region to the bottom right. We can now ver- 
ify our tentative identifications of these stable regions 
of the diagram. We can also determine that the unsta- 
ble fixed point controls the transition between the glass 
and the superfluid, the superfluid-insulator transition of 
our disordered rotor model. This allows us to draw the 
schematic picture of the flow diagram that we presented 
in Figure 1. 

We will now focus on the critical region and extract 
critical exponents governing the glass-superfluid transi- 
tion. Estimating these exponents requires formulating 
scaling hypotheses for the behavior of physical quantities 
in the critical region. In the case of s n 



28. 



(70) 



s m ax = I^i (jj =L d f~h 1 (Li(g-g c )) 

Exactly at criticality, <s max ~ L df asymptotically, so plot- 
ting vs. (g — g c ) generates a crossing of the curves for 
different lattice sizes. Slightly away from criticality, the 
power law behavior should persist if L < £. For L > £, 
the system recognizes that it is not critical and we should 
see a crossover to logarithmic growth on the glassy side 
and L 2 growth on the superfluid side. Hence, j is the im- 
portant ratio near criticality, and this motivates scaling 
form (70). The scaling form, in turn, implies that we can 
produce scaling collapse by plotting vs. (g — g c )L» . 
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This is what we have done in panel (a) of Figure 16. This 
scaling collapse leads to the estimates: 

i/« 1.09 ±0.04 (71) 

d/«1.3±0.2 (72) 

These estimates and error bars are obtained through 
the following procedure: first, to find an estimate of the 
critical value of the tuning parameter g, we examine the 
behavior of the sample average of 52 vs. g. Since 52 varies 
as a power law in L exactly at criticality and grows more 
slowly within the phases, S2 should exhibit a maximum 
at g c . We can find the error Ag c on our estimate of g c 
by taking the window of values of g around g c for which 
the sample average of 52 is within one error bar of the 
maximum. To proceed to obtain estimates for v and df, 
we now partition our data into bins of size samples. 
For example, the first bin may consist of the first = 50 
samples for each value of the tuning parameter and each 
lattice size L. Our immediate goal is to find the best 
values of v and df for this subset of our data. We first 
focus on the data for g — g c and do linear regression 
to find the best exponent that describes the power law 
growth of the average value of <s max with L. This gives 
an estimate of df for the bin. Next, for two lattice sizes 
(typically, L = 150 and L — 300), we compute an average 
of 5 max over the samples in the bin for several values of 
the tuning parameter around g c . Then, using the previ- 
ous estimate of df for the bin and the scaling hypothesis 

(70) , we test various values of v until we achieve the best 
collapse. Thus, we also obtain an estimate of v for the 
bin. From the distribution of estimates for the different 
bins, we can find mean values and error bars for df and 
v. However, these error bars do not take into account the 
error on our estimate of the critical point. To propagate 
this error, we need to repeat the critical point estimation 
procedure using g c + Ag c and g c — Ag c as our estimates 
of the critical point. 

We have repeated the scaling collapse of 5 max for a 
different choice of initial distributions: bimodal P%{U) 
and Pi (J) oc J -1,6 . Ultimately, we have been able to 
recover estimates for v and df which are consistent with 

(71) and (72): 

i/« 1.1 ±0.1 (73) 

df& 1.2 ±0.2 (74) 

A completely analogous scaling hypothesis can be 
made for 52: 

s 2 =L d fh 2 (jj =L d fh 2 (Li(g-g c )) (75) 

Then, the exponents needed to produce the collapse in 
panel (b) of Figure 16 are: 

i/« 1.06 ±0.09 (76) 



d/« 1.31 ±0.07 (77) 

Since all the estimates (71)-(77) are consistent, we will 
proceed using our tightest estimates on these exponents: 
(71) and (77). 

We now note that a scaling hypothesis can also be 
formulated for the density p(s). From equation (45), we 
see that, for fixed lattice size L, s T p(s) should depend 
only on the combination | ~ s(g — g c )°- Hence, by 
plotting s T p(s) vs. s(g — g c )i and tuning until the curves 
for different choices of g collapse, we ought to be able to 
extract estimates for r and a. On the other hand, r and 
a are related to v and df through the scaling relations 
(61) and (57), so from the estimates (71) and (77), we 
can infer: 

r ^ 2.53 ±0.08 (78) 

a « 0.70 ± 0.04 (79) 

In Figure 20, we attempt to produce collapse of p(s) for 
L = 300 lattices using these estimates of r and a. To 
produce this plot, we need to discard data points for small 
cluster sizes (s < 30), because they are non-universal, 
and for large cluster sizes (s > 100), because they are 
noisy. Once we do this, the collapse works fairly well, 
indicating that we have found a consistent set of critical 
exponents obeying the necessary scaling relations. 

Comments on the Insulator- Insulator Transition 

Before proceeding further, we should note that our nu- 
merics do not accurately capture the boundary between 
the Mott insulator and the Mott glass. Several authors 
have argued that we should expect glassy behavior to 
occur whenever the ratio of the largest bare Josephson 
coupling to the lowest bare charging energy exceeds 
the value of the clean transition, because this condition 
allows for the presence of regions in which the system 
locally looks like it is in the superfluid phase 6 ' 15 . How- 
ever, in the strong disorder RG treatment, some distri- 
butions that satisfy this criterion still produce Mott in- 
sulating behavior out to the largest lattice sizes that we 
investigate. Since the glassy phase occurs due to rare- 
regions or Griffiths effects, a finite size system will only 
look glassy if it is large enough for the rare-regions to 
appear. This suggests that some choices of parameters 
which yield Mott insulating behavior on finite size lattices 
may actually yield glassy behavior in the thermodynamic 
limit. Of course, this difficulty necessarily afflicts all nu- 
merical methods (except those that rely on mean-field 
type approximations 29 ), since they are confined to oper- 
ate on finite size systems. 

We will not comment on this transition further, but 
we take this opportunity to refer the reader to papers by 
Kriiger et al. and Pollet et al., which present two view- 
points on the Mott insulator to Bose glass transition 6 ' 30 . 
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C. Identifying the Glass 

Finally, we return to the question of the identity of 
the glassy phase. Is the phase a Bose glass or a Mott 
glass? A definitive diagnosis requires a measurement of 
the compressibility: 
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dp 



(80) 
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The compressibility is more subtle to measure than the 
quantities that we have already discussed. Any finite 
size system is gapped and therefore incompressible. On 
the other hand, in the thermodynamic limit, the gap can 
vanish and the compressibility need not be zero. 

How can we measure the compressibility of the glassy 
phase in the RG? In Figures 17 and 18, we presented 
data for the density (number per unit area) of clusters 
with a given gap A. With this density profile in hand, we 
can calculate the density of particles introduced to the 
system by a small chemical potential shift p: 
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Here, n(A) is the number of particles added to a clus- 
ter with gap A if the chemical potential is p. If p(A) 
stays finite as A — >> 0, the integral is divergent, and the 
system is infinitely compressible. Suppose alternatively 
that p(A) vanishes as A^ for small A. Then: 
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The derivative of the integral vanishes at p = for j3 > 
0. Thus, the system is incompressible. Comparing to 
the data shown in Figures 17, we see that there is no 
numerical evidence for a finite glass compressibility; the 
gap density appears to vanish even faster than a power 
law as A —> 0. This is consistent with the behavior of 
p(s) in equation (45), since A is expected to scale as s" 1 . 
Hence, the numerics imply that the Mott glass intervenes 
between the Mott insulator and the superfluid in this 
model. 

At first glance, the preceding argument may be trou- 
bling. Due to the shrinking of the low A peak in panel 
(b) of Figure 18, the gap density p(A) also appears to 
vanish as A —> in the superfluid phase. The caveat is 
that it is necessary to more carefully evaluate the com- 
peting effects of the shrinking and the propagation. The 
low A peak in Figure 18 represents the macroscopic su- 
perfluid clusters that form in the superfluid phase. These 
clusters do not appear in proportion to the surface area 



of the sample, as is the case for rare-regions clusters; in- 
stead, one such macroscopic cluster appears in each sam- 
ple. Therefore, the density of macroscopic clusters will 
go as jp. , and this is responsible for the shrinking of the 
low A peak. The propagation of the peak, meanwhile, 
reflects the fact that the gap closes as L~ 2 . For a fixed 
choice of /i, the number of bosons that will be added to 
these macroscopic clusters scales as: 
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for large L. Then, the density of particles introduced to 
the system is: 



(84) 



This directly implies that the compressibility (80) is a 
constant in the thermodynamic limit, so we recover the 
expected result that the superfluid phase is compressible. 



VI. CONCLUSION 

While the interplay of disorder and interactions in 
bosonic systems has attracted considerable interest for 
nearly three decades, the random boson problem remains 
a fertile source of intriguing physics. In this paper, 
we have investigated a particular model of disordered 
bosons, the two-dimensional rotor (or Josephson junc- 
tion) model. Our strong disorder RG analysis suggests 
the presence of an unstable finite disorder fixed point 
of the RG flow, near which the coupling distributions 
flow to universal forms. Furthermore, the strong disorder 
renormalization procedure indicates the presence of three 
phases of the model: the Mott insulating and superfluid 
phases of the clean model are separated in the phase dia- 
gram by an intervening glassy phase. The unstable fixed 
point governs the transition between the superfluid and 
this glassy phase, and the transition is driven by a kind 
of percolation. The RG procedure also provides evidence 
that this glassy phase is, in fact, the incompressible Mott 
glass. 

Our work is a numerical extension into two dimen- 
sions of the one-dimensional study by Altman, Kafri, 
Polkovnikov, and Refael 13 . The 2D fixed point, how- 
ever, differs from the ID fixed point in an important 
way. The ID fixed point occurs at vanishing interac- 
tion strength (all charging energies Uj = 0). Thus, it 
corresponds to a completely classical model and reveals 
that the superfluid-insulator transition can be tuned by 
varying the width of the Josephson coupling distribution 
at arbitrarily small interaction strength. The 2D fixed 
point is, in contrast, fully quantum. Indeed, in the crit- 
ical distributions generated by the strong disorder RG, 
the charging energy distribution is peaked near the RG 
scale while the Josephson coupling distribution is peaked 
well below. 
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On the other hand, the fixed point that we have iden- 
tified in this paper is similar to its ID counterpart in 
that it does not exhibit infinite randomness. Finite dis- 
order fixed points are not optimal settings for strong dis- 
order renormalization analyses, because the procedure 
does not become asymptotically exact near criticality 
and is, in this sense, an uncontrolled approximation. We 
have proceeded with such an analysis anyway. In doing 
so, we have found a robust fixed point controlling the 
superfluid-insulator transition and phases exhibiting rea- 
sonable physical properties. While this may be surpris- 
ing given the perils of the method, we have attempted 
to argue for the appropriateness of the method, as an 
approximation, through an analysis of the RG steps in 
light of the forms of the fixed point distributions P U nw(U) 
and Puniv(^)- We certainly acknowledge that there are 
other subtleties due to the lack of infinite randomness; 
for example, the notion of a superfluid cluster is not 
completely sharp, and consequently, percolation of su- 
perfluid clusters can only be an approximate picture of 
the transition 31 . Nevertheless, the structure of the fixed 
point Josephson distribution (36) suggests that the pic- 
ture may be a good approximation, and we take this 
opportunity to remind the reader that we extensively 
discuss the reliability of the RG, in light of the prop- 
erties of the fixed point, in Appendix D. Moreover, the 
self-consistency of our numerical results, especially the 
striking universality and robustness of the unstable fixed 
point, gives us confidence that our strong disorder RG 
analysis provides useful information about the system. 
With the potential caveats in mind, we therefore turn to 
exploring connections with other theoretical, numerical, 
and experimental work. 

The Mott glass phase of the two-dimensional model is 
the straightforward analogue of the phase found in one 
dimension by Altman et al. The charging gap, the energy 
needed to add or remove a boson from the system, van- 
ishes due to the presence of arbitrarily large rare-regions 
of superfluid ordering. However, there is no true long 
range order because these rare-regions grow subexten- 
sively with system size. If a small chemical potential 
shift is turned on, it becomes energetically preferable to 
add bosons somewhere in the system, specifically in the 
largest of the rare-region superfluid clusters. Neverthe- 
less, these clusters do not occur with sufficient number to 
produce a finite density of bosons, and the glass remains 
incompressible. In a Monte Carlo study of a model sim- 
ilar to ours, Prokof'ev and Svistunov previously found 
evidence for a glassy phase in which the compressibility 
vanishes for this reason 32 . Moreover, the Mott glass that 
was identified by Roscilde and Haas in a related spin-one 
model also relies on the same mechanism 33 . The origi- 
nal proposal of Giamarchi, Le Doussal, and Orignac is, 
however, fundamentally distinct 12 . In their Mott glass, 
the charging gap remains finite, guaranteeing a vanishing 
of the compressibility; however, gaplessness is achieved 
through the closing of a mobility gap for transport of par- 
ticles between nearby insulating and superfluid patches. 



Sengupta and Haas have argued that particle-hole sym- 
metry, a crucial ingredient in the formation of our Mott 
glass, is not necessary for the realization of the phase 
through this alternative mechanism 34 . 

In the superfluid phase, true long range order emerges 
because the largest cluster scales with the size of the lat- 
tice. In this sense, this cluster is macroscopic. Despite 
this, near the transition, the macroscopic cluster may 
only occupy a small fraction of the total number of lat- 
tice sites. Because the clustering procedure can merge 
sites that are not nearest neighbors in the bare lattice, 
the fraction of insulating sites may actually exceed the 
standard 2D square lattice percolation threshold. Even 
with such a large fraction of insulating sites, a superfluid 
phase can still exist because virtual tunneling processes 
can carry super current through the insulating sites, al- 
lowing for macroscopic superfluidity on the "depleted" 
lattice that forms when the insulating sites are removed 
from the lattice by site decimation. 

Nevertheless, the Mott glass to superfluid transition 
of our model should be contrasted with transitions that 
arise when disorder is introduced to a 2D square lattice 
by bond or site depletion 33 ' 35 39 . Models of the latter 
type only have the opportunity to form long range or- 
dered phases when the underlying lattice is percolating. 
This percolation is purely classical and exhibits all the 
critical properties expected of standard site or bond per- 
colation on a square lattice 28 . The superfluid-insulator 
transition is, in general, distinct from this transition; 
once the underlying lattice percolates, the bosonic model 
defined on that lattice may still exhibit Mott insulating, 
glassy, and superfluid phases 33 . In contrast, the only 
percolation-type process in our model is the one that ac- 
tually drives the superfluid-insulator transition. The crit- 
ical properties of this transition differ dramatically from 
those of classical 2D square lattice percolation, because 
the transition is not a purely geometric phenomenon. In- 
stead, there are quantum tunneling processes overlaid on 
top of a geometric structure 40 . 

We have remarked that several experimental groups 
are currently working on studying disordered bosonic sys- 
tems in ultracold atoms 8 11 . Dirty and granular super- 
conductors provide another experimental context that 
may be relevant to the physics of this paper. The 
question of the applicability of bosonic pictures to the 
2D superconductor-insulator quantum phase transition 
is long-standing. Recently, this problem has been ad- 
dressed numerically by Bouadim et al., who used quan- 
tum Monte Carlo simulations to study a fermionic model 
of the transition and found that bosonic physics emerges 
at criticality 41 . This issue has also been addressed exper- 
imentally by Crane et al. These authors studied indium 
oxide thin films and measured a superfluid stiffness in 
the insulating state, indicating that Cooper pairs may 
survive into the insulating region 42 . This leaves open 
the possibility that the transition is driven by percola- 
tion of superconducting regions, a possibility that has 
also been explored in experiments on granular super- 



24 



conductors by Frydman et al. and Sherman et al. 43 ' 44 . 
During the preparation of this manuscript, we learned 
of a recent experiment by Allain et al. on tin-decorated 
graphene. Intriguingly, this experiment may point to a 
superconductor-insulator transition that is bosonic in na- 
ture, driven by percolation, and characterized by critical 
exponents similar to those identified by our work 45 . 

Motivated by the link between the Mott glass and 
particle-hole symmetry, Roscilde and Haas have proposed 
a different class of experimental systems in which Mott 
glass physics may be present: nickel based spin-one an- 
tiferromagnets. The advantage in these spin systems is 
that it may be easier to realize particle-hole symmetry in 
the guise of Z 2 symmetry 33 . Very recently, Yu et al. have 
followed up on this proposal with an experimental inves- 
tigation of Bose and Mott glass phases in bromine-doped 
dichloro-tetrakis-thiourea-nickel 46 . This system is three- 
dimensional, so the character of the transition would nat- 
urally differ from what we have calculated in this paper. 

One immediate extension of our study would be to con- 
sider adding random filling offsets to the disordered rotor 
model, as Altman et al. did in a follow-up to their work 
on the ID model 14 . The intuition from ID suggests that 
such a model would contain a Bose glass phase. On the 
other hand, very recent Monte Carlo work by Wang et 
al. indicates that the Mott glass survives the substitution 
of exact particle-hole symmetry with statistical particle- 
hole symmetry 47 . In one dimension, Altman et al. found 
that the universality class of the transition (but not the 
identity of the glassy phase) is independent of the sym- 
metry properties of the random filling offsets, but the 
situation may differ in d > 1; this remains to be under- 
stood. 

Another interesting extension may be to study the ro- 
tor model defined on random networks. Suppose we do 
not begin with a square lattice but rather with a gener- 
alized network of mean connectivity z = 4. At its critical 
point, would such a model flow to the same fixed point as 
the model defined on a square lattice? The fact that the 
strong disorder RG modifies the initial lattice structure 
into a more general network suggests that this may be the 
case for at least some types of random networks. Next, 
suppose we vary the mean connectivity from z = 4. Is 
there a range of connectivities for which random network 
models access our fixed point? 

Perhaps it would be better to precede such investiga- 
tions with a better characterization of the fixed point 
itself. In one dimension, Altman et al. were able to write 
down master equations for the RG flow, solve them to 
find fixed point distributions, and then verify numerically 
that these distributions are stable 13 . In two dimensions, 
a direct analytical approach is more difficult, and it re- 
mains to be seen whether such an approach is tractable. 
Our work provides suggestive numerical evidence regard- 
ing the forms of the universal distributions that char- 
acterize the critical point of the disordered rotor model 
(7). A Lyapunov analysis of these distributions, in which 
the RG is used as a tool to identify irrelevant directions 



in the space of possible distributions, could be a useful 
step in clarifying the critical forms still further. Then, 
analytically verifying these forms as attractor solutions 
of the RG flow may be an easier task than analytically 
identifying them would have been in the absence of any 
numerical guidance. 

Recent work by Vosk and Altman suggests yet another 
direction for connecting the results of the RG to exper- 
iment. These authors have derived the ID version of 
the rotor model as an effective description of continuum 
bosons. In doing so, they have established a connection 
between the strong disorder RG treatment of Altman, 
Kafri, Polkovnikov, and Refael and cold atom experi- 
ments on rubidium-87. Remarkably, the distributions 
that Vosk and Altman derive are of the same form as 
the fixed point distributions found by the strong disor- 
der RG 13,20 . If such a treatment can be extended to the 
2D case, that would be very valuable, both as a clarifica- 
tion of the critical behavior and as an indication of the 
relevance of the results to current experiments. 
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Appendix A: Sum Rule vs. Maximum Rule 

In this appendix, we present a short argument for why 
it may be preferable to use the sum rule (see equations 
(32) and (33)) instead of the maximum rule (see equa- 
tion (34)). In dimensions greater than one, it should be 
easier to form ordered (e. g. superfluid) phases. Indeed, 
the transition for the clean rotor model occurs when jj is 
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substantially smaller than one 22 . Suppose we begin with 
the clean model at its critical point and then disorder it 
by increasing some Josephson couplings and decreasing 
some charging energies. Suppose further that we do this 
such that that the greatest Josephson coupling J max is 
still less than the weakest charging energy U m { n . Then, 
using the maximum rule, the strong disorder renormal- 
ization procedure will predict no cluster formation at all. 
In other words, it will predict the ground state to be a 
Mott insulator, and this result seems inconsistent with 
the location of the clean transition. As mentioned previ- 
ously, several authors have argued that, if exceeds 

the ratio jj at the clean transition, then we should ex- 
pect glassy behavior, because there can be rare-regions in 
which the system looks locally superfluid 6 . With the sum 
rule, the procedure has a mechanism to circumvent this 
inconsistency. The Josephson coupling scale can actually 
grow through decimation and compete with the charging 
energy scale. Thus, there can be cluster formation, and 
the procedure can find ground states that are glassy or 
superfluid, even when all Josephson couplings of the bare 
model are less than the minimum bare charging energy. 
This indeed occurs in the numerics, as we have noted 
while presenting the numerical data above. 

The notion of the Josephson coupling scale increasing 
through the RG may be a source of concern to some read- 
ers. The increase actually corresponds to the generation 
of multiple effective couplings between two sites through 
different paths in the lattice. This still happens when 
the maximum rule is used, but it is hidden through the 
discarding of certain paths. If the coupling through each 
path is treated as an independent Josephson coupling, 
then the Josephson coupling scale does decrease as the 
renormalization proceeds. However, when it is time to 
determine the next decimation step, it is necessary to 
consider all of the couplings between any two sites. For 
this reason, it makes sense to sum all the couplings into 
one effective coupling between the sites. 



Appendix B: Measuring Physical Properties in the 

RG 



Here, we work out two examples of how estimates of 
physical properties can be extracted from the RG proce- 
dure. 



Particle Number Variance 



First, consider the quantity: 



(Bl) 



This particle number variance gives the mean squared 
number fluctuation away from the large filling, summed 
over all sites in the lattice. When normalized by , we 



find numerically that this quantity stays constant as the 
system size is increased for all choices of distributions 
and parameters. As such, this quantity is completely 
uninteresting for discriminating between phases of the 
model, but we do calculate it for comparison to exact 
diagonalization in Appendix D. 

The calculation of the particle number variance (Bl) 
is most straightforward when clusters do not form, so 
let us first consider the case where some site X is not 
clustered with any other site during the renormalization. 
At some stage in the procedure, the site is decimated 
away. The number fluctuation on site X is locked to 
zero at leading order with corrections incorporated in 
second order perturbation theory. An approximation to 
the ground state value of (n 2 x ) can be obtained from the 
perturbative expansion of the state. In particular: 
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When clusters do form, the calculation is trickier, but 
it can be performed by carefully keeping track of how the 
operator that we are targeting is written in terms of the 
cluster and relative number operators introduced in link 
decimation. To illustrate this, suppose we are trying to 
measure the operator: 
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The factors aj and a k are just numbers. In (Bl), all 
these factors are one, but we will motivate the inclusion 
of more general a factors here shortly. If sites j and k are 
merged into a cluster, then we switch to the operators nc 
and ur to find: 
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During link decimation, the relative coordinate is spec- 
ified, so the expectation value of the final term can be 
found immediately from the harmonic approximation: 



(dj +a k )(h 2 R ) « (dj +a fc ) 
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(B5) 



Furthermore, the harmonic theory also predicts that the 
expectation value of the term linear in fiR will vanish. 
The calculation of the term proportional to r? c must be 
deferred to later in the renormalization procedure. Thus, 
we keep the operator h 2 ^ in the portion of the sum (Bl) 
that remains to be evaluated, where it appears just like 
the ft 2 for a bare site, but multiplied by a renormalized 
ac coefficient. This was the motivation for including the 
a factors; though the bare values of these factors are all 
equal, different values can be generated through cluster 
formation. If the cluster is merged with another clus- 
ter in a future link decimation, we repeat the procedure 
above. When the cluster is finally decimated in a site 
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decimation, the cluster's contributions to the sum (Bl) 
are calculated through equation (B2) and then multiplied 
by the appropriate a factor. 

The procedure outlined above can run into a difficulty 
that we can anticipate by thinking about the two-site 
problem. Suppose two sites, labelled 1 and 2, are con- 
nected by a Josephson coupling J12. Furthermore, sup- 
pose Ui > U2, but both charging energies are greater 
than J 12 . Then, we would decimate site 1 first and ob- 
tain an estimate for the site's contribution to the particle 
number variance (Bl) through equation (B2). Next, we 
would decimate site 2 and find that it does not contribute 
at all to (Bl) because there are no remaining links. How- 
ever, the contribution of site 2 should, in fact, be equal 
to that of site 1, so our estimate is off by a factor of 
two. We can verify this by adopting cluster and relative 
coordinates (23) and then calculating (Bl) by doing per- 
turbation theory for the relative coordinate. To partially 
resolve this difficulty, we can keep track of which sites are 
unclustered by the RG process. At the end, we can return 
to the original lattice and calculate the contributions of 
these unclustered sites to (Bl) using the bare couplings. 
For sites that are clustered by the RG, we reason that 
the main contribution to (Bl) comes from internal fluc- 
tuations of the cluster, so this correction may not be so 
important. 

It is possible to raise another objection to our calcula- 
tion of (Bl). The perturbation theory that leads to the 
result in (B2) incorporates the charging energies on sites 
neighboring site X. However, the perturbation theory 
leading to the RG rule (21) does not. How do we resolve 
this contradiction? When we calculate effective Joseph- 
son couplings in the RG, what we want to calculate is 
the effective rate of tunneling, once a boson has left one 
neighboring site and before it arrives at another, through 
the link-site-link system. For this purpose, it is appro- 
priate to treat the sites neighboring site X as fictitious 
charging energy-free islands. On the other hand, when 
calculating observable quantities like (Bl), it is impor- 
tant to account for the fact that the ability to move a 
particle or hole from site X to a neighboring site also 
depends on how hard it is to charge the neighboring site. 

2. Susceptibility 

In our studies of physical properties of the phases of the 
disordered rotor model, we considered the susceptibility 
X- We defined this quantity in equation (41) as the linear 
response of the system to a uniform ordering field (40) 
coupling to cos(0j). 

To calculate x> we consider how to calculate the ex- 
pectation value: 

5><cos(^)> (B6) 

3 

in the presence of an infinitesimal ordering field h. As 
with the a-factors in the calculation of the particle num- 



ber variance, all the bare bj — 1. When clusters form, 
the effective 6-factor for the cluster can differ from unity. 

If a bare site is decimated, then perturbation theory in 
h gives: 

(cos(^)) = A (B7) 

Since bx = 1 for a bare site, (B7) is the contribution of 
site X to the quantity (B6). 

When a link decimation joins two sites into a cluster, 
the corresponding terms in the sum (B6) merge as: 

bj cos(4>j) + b k cos(0/e) = bj cos(4> c + H&r) 

+b k cos(^c - Vk^R) 

« bjCos(4>c)(cos(fij4> R )} (B8) 
-bj sm(4) C )(sin(fij4) R )) 
+b k cos(0 c )(cos(/i /e jR )) 
+b k sin(0c)(sin(/i fe 0i?)) 

= (bjCnwj + b k c DW ,k) cos((/>c) 

In this calculation, \ij and \i k are the ratios introduced 
in equation (30), and cdwj and cnw,k are precisely the 
Debye- Waller factors given in equation (31). We can read 
off the renormalized 6-factor for the cluster from the cal- 
culation above, and the resulting expression is given in 
equation (42). 

Next, it is important to note that the ordering field 
terms in the Hamiltonian (40) transform in the same 
way. In other words, the term coupling to cos (0c) in 
the Hamiltonian should appear multiplied by be after 
a merger. Physically, this corresponds to the fact that, 
when the cluster phase rotates, sc bare phases rotate 
semi- coherently. Complete coherence is lost due to quan- 
tum fluctuations, which are accounted for by Debye- 
Waller factors. The factor be can be thought of as the 
effective number of rotors that are coherently ordering 
with the field. The energetic cost of <\>c straying from 
the direction of the ordering field is therefore amplified 
by this factor. When the cluster is finally decimated, the 
perturbation is amplified by this amount. Furthermore, 
cos(^c) appears in the sum (B6) multiplied by be, so the 
total contribution of the cluster to the sum is: 

hh 2 

bc(cos^c)) = (B9) 
Uc 

From equations (B7) and (B9), the calculation of the 
linear response to an infinitesimal h (i.e. the susceptibility 
x) follows immediately. 

Appendix C: The Pi(U) Gaussian, Pi(J) Power Law 
Data Set 

In the work reported above, we have explored the 
superfluid-insulator transition of the disordered rotor 
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model using several different species of disorder (see 
equations (9)-(12)). In doing so, we have exposed uni- 
versal features of the critical behavior. After providing 
numerical evidence of this universality however, we focus 
on data for one particular choice of initial distributions. 
In this appendix, we describe this choice of distributions 
in greater detail. 

The choice of distributions in question first appears in 
Figure 7. The initial distribution of charging energies 
Pi(U) is taken to be a Gaussian with center at Uo and 
width gjj = 2. Hence, the form of the distribution is: 

t/o) 21 



P%(U) oc exp 



(U 



(CI) 



Recall that the Gaussian distribution is truncated at 
3<Ju = 6, so the distribution only has weight in the in- 
terval U G (Uo — 6, Uo + 6). We leave Uo unspecified for 
the moment, because it is the parameter that we use to 
tune through the transition. 

The initial Josephson coupling distribution P%(J) is a 
power law of the form J -16 . We choose the cutoffs so 
that the distribution has weight for J G (0.5, 100). This is 
a very wide power law distribution, so the corresponding 
flows begin well above the unstable finite disorder fixed 
point in the numerical flow diagrams. Explicitly: 

Pi(J) «0.413J -1 ' 6 (C2) 

For this choice of distributions, we have acquired data 
for Uo = 400, 20, 18, 16, 14, 12, 10, 9.8, 9.6, 9.4, 9.2, 8.8, 
8.6, 8.4, 8.2, 8, 7.5, 7, and 6.5. To more finely target the 
transition, we have probed the interval 9.1 > Uo > 8.9 
in increments of 0.01. We have always acquired data for 
L = 25, 50, 75, 100, 150, 200, and 300. In all cases, we 
have pooled data for 10 3 disorder samples. The peak in 
the data for S2 vs. Uo gives the following estimate of the 
critical point: 



Uo, c -8.97 ±0.02 



(C3) 



Close to criticality, it is better to use a lower value of the 
thresholding parameter a. For several values of Uo such 
that 10 > L^o > 8, we have run the RG with = 10 -5 
and o^ = 5xl0 -6 to test for convergence. Further away 
from criticality, we have instead used ^ = 5x 10 -5 and 
an — 2.5 x 10 -5 . Figure 22 shows a test of the convergence 
of the maximum cluster size s max in the thresholding 
parameter a. We plot: 

°max 



(C4) 



vs. the tuning parameter Uo- The ratio (C4) is essentially 
always within two error bars of unity. We take this as an 
indication that physical properties have converged. 



Appendix D: Further Comments on the Use of the 
Strong Disorder RG in a Finite Disorder Context 

This appendix is devoted to exploring, in further de- 
tail, the validity of the RG procedure. We first expand 




FIG. 22: A test of the convergence of s ma x in the thresholding 
parameter a. The variable v is the ratio of the estimate of 
Smax for an, the less conservative value of the thresholding 
parameter, to ct£, the more conservative value. We plot v 
against the tuning parameter Uo. The closest data point to 
criticality (Uo = 8.97) is indicated with a diamond. Note that 
v > 1 typically indicates convergence since less conservative 
thresholding (higher a) corresponds to throwing away more 
bonds and, therefore, biases the system away from cluster 
formation. Smaller values of an and ae are used in the vicinity 
of the transition. See the text of Appendix C for details. 



upon the argument, introduced earlier in the paper, for 
the reliability of the RG near criticality. Then, we move 
away from criticality and assess the performance of the 
RG in the various phases of the disordered rotor model 
(7). Next, we focus on each of the RG steps, consider 
circumstances in which they may fail to capture impor- 
tant physics, and formulate tests to ensure that the RG 
is trustworthy in these situations. Finally, we present a 
comparison of the RG to exact diagonalization of small 
systems. 



1. Review of the Argument for the RG at 
Criticality 

Our confidence in the RG procedure near criticality 
rests on the form of the critical Josephson coupling dis- 
tribution, reported in equations (36) and (39). Infinite 
randomness develops when P(J) oc J -1 , and our nu- 
merical evidence suggests that the critical distribution of 
the disordered rotor model decays even more strongly 19 . 
Nevertheless, the critical distribution does not exhibit 
infinite randomness, because as seen in Figure 11, it is 
cut off from below. Recall that the lower cutoff of the 
"Josephson coupling distribution" is set by our choice 
to retain, in statistics, only the dominant 2N Joseph- 
son couplings, where N is the number of sites remain- 
ing in the effective lattice. Then, the appropriate way 
to interpret the distributions in Figure 10 is the follow- 
ing: penetrating any given site in the effective lattice, 
there are likely to be on the order of four Josephson cou- 
plings drawn from the depicted distribution. There may 
be other links penetrating the site, but these will be even 
weaker. We can estimate the typical strength of the four 
dominant links by comparing the median of the critical 
Josephson coupling distribution to the maximum. For 
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the closest approach to criticality with the initial distri- 
butions described in Appendix C, we find that the ratio 
is approximately 0.11 ± 0.01 near the fixed point. 
Hence, the typical link is quantitatively weak compared 
to jj ~ 0.345 at the clean transition 22 . 

The considerations above form a strong argument for 
the validity of the site decimation RG step. Here, we 
seek out the dominant effective charging energy in the 
effective lattice, and treat the links penetrating the site 
in perturbation theory. This perturbation theory is likely 
a very good approximation, because the Josephson cou- 
plings penetrating the site in question are usually ex- 
tremely weak. 

Now, we turn to the link decimation step, in which we 
seek out the largest Josephson coupling in the lattice and 
merge the corresponding sites into a cluster. Although, 
the Josephson coupling being decimated is the largest en- 
ergy scale in the system, there is a high probability that 
all the other links penetrating the two sites being joined 
will be very small. However, the critical distribution of 
charging energies is not peaked at low U. The structure 
of the distributions plotted in Figure 12 suggests that 
it is quite likely that one or both sites being merged will 
have a charging energy of the same order as the RG scale, 
thus violating the strong disorder hypothesis. Our treat- 
ment of the quantum fluctuations of the relative phase of 
the cluster is based on the harmonic approximation (27). 
Is this approximation appropriate when the charging en- 
ergies of the two-site problem are comparable in magni- 
tude to the Josephson coupling? Alternatively, do the 
quantum fluctuations grow so large that the clustering 
becomes meaningless? We address this question as fol- 
lows: using the fact that the remaining links are weak, we 
isolate the two site problem and solve it exactly, treating 
the remaining links via second order perturbation the- 
ory. Comparing the results of the RG with the exact 
solution, we find that, even in this worst case scenario, 
the RG produces reasonably accurate effective couplings. 
The evidence for this claim is given in section 3 below. 



2. Reliability of the RG in the Phases 

How reliable is the RG when we move away from criti- 
cality and into the phases of the disordered rotor model? 
To gain some insight into this question, we can consult 
Figure 23, which expands upon Figures 11 and 13 by 
plotting renormalized J and U distributions away from 
criticality. 

Proceeding into the glassy regime, the arguments pre- 
sented above for the validity of the RG near criticality 
generally become stronger. The primary reason for this 
is that the renormalized J distributions become progres- 
sively broader than at criticality. In the flow diagrams 
(Figures 6-9), this is reflected in the apparent divergence 
of ~. Consequently, the assumption of isolating local 
degrees of freedom becomes better as we get deeper into 
the glass. One complication is that the renormalized U 




ln(U/Q) 

FIG. 23: In panel (a), a sweep of renormalized J distribution 
through the glass and into the superfluid. The initial distribu- 
tions are those described in Appendix C: P%(U) is Gaussian 
and Pi(J) oc J -1 " 6 . All data is taken for L = 300 lattices, 
and the renormalized distribution is computed when 300 sites 
remain in the effective lattice. The values of Uo shown are 
18, 12, 9.6, 8.97 (near-critical), 8.4, and 7.5. Panel (b) shows 
the corresponding sweep of the renormalized U distribution 
at the same stage of the RG. 



distribution becomes more strongly peaked near the RG 
scale. This may pose trouble for the link decimation step 
and makes it especially important to consider the reli- 
ability of this step when a strong Jjk couples two sites 
with charging energies Uj ~ Uk ~ Jjk- As mentioned 
previously, we will study this worst-case scenario in part 
3 of this appendix and find that the RG still works rea- 
sonably well. 

Now, we turn to the superfluid phase. As we proceed 
away from criticality, the renormalized charging energy 
distribution becomes flatter and broader. The broad- 
ening of this distribution implies that the likelihood of 
encountering a strong Josephson coupling that connects 
sites with comparably strong charging energies decreases 
as we get deeper into the superfluid phase. However, 
the J distribution also seems to become flat deep in the 
superfluid, and this is problematic. For example, dur- 
ing link decimation, it may not be a reasonable approx- 
imation to isolate the two-site problem centered on the 
strongest Josephson coupling. As we proceed into the 
superfluid, it is necessary to be more dubious of the RG 
results; nearer to criticality however, the arguments used 
at the critical point are probably approximately valid. 

The plots in Figure 23 attempt to elucidate system- 
atics in the behavior of the renormalized distributions 
in the insulating and superfluid regimes, but we should 
mention that this figure should be interpreted with some 
care. For the choice of distributions in Appendix C, flows 



terminating in the insulating or superfluid regions of the 
flow diagram nevertheless propagate in the vicinity of the 
unstable fixed point along the way. This can be seen, for 
instance, in Figure 7. For certain choices of initial dis- 
tributions, there can be flows towards the insulating or 
superfluid regimes which never propagate anywhere near 
the unstable fixed point. Consequently, the RG never 
has an opportunity to wash away the details of the ini- 
tial distributions and allow the universal properties of the 
fixed point to emerge. Hence, the renormalized distribu- 
tions generated along such flow trajectories are unlikely 
to exhibit the clean systematic properties seen in Figure 
23. 
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FIG. 24: The graph structure for the tests reported in Figures 
25 and 26. The links J13, J14, and J25 are assumed to be 
perturbatively weak. The charging energies U\ and U2 and 
the Josephson coupling J 12 are varied to explore potentially 
troublesome scenarios. 



3. Analysis of Potential Problems with the RG 

Next, we address some potential difficulties with the 
arguments presented above for the reliability of the renor- 
malization procedure in the critical region. These diffi- 
culties are rooted in the lack of strong randomness in the 
critical charging energy distribution. 

Consider the lattice geometry shown in Figure 24. Sup- 
pose that all links except the Josephson coupling between 
sites 1 and 2 are perturbatively weak. Now, suppose fur- 
ther that J12, E/i, and U2 are comparable in magnitude, 
but J 12 is the largest of the three. We would like to for- 
mulate a test of whether the RG appropriately handles 
this situation. Within the RG, a link decimation would 
merge sites 1 and 2 into a cluster. All links penetrat- 
ing sites 1 and 2 would be modified by their correspond- 
ing Debye- Waller factors (31) cdw,i an d cdw,2- Then, 
because all the remaining links are assumed to be very 
weak, the cluster of sites 1 and 2 will be decimated, pro- 
ducing an effective coupling between sites 3 and 5: 



J35,RG = CDW,lCDW,2 



(Dl) 



Another approach to calculating this effective coupling 
would be the following: take the two-site problem of sites 
1 and 2 and exactly diagonalize it. Then, to leading 
order, sites 1 and 2 should be locked into the two-site 
ground state, with perturbative corrections coming from 
the Josephson couplings J13, J14, and J25. This pertur- 
bation theory leads to an effective coupling through the 
sites 1 and 2. This alternative procedure is perhaps more 
appropriate, because it does not presuppose the harmonic 
approximation. In Figure 25, we assess how much of an 
error we make by adopting the harmonic approximation. 
Holding J fixed, we sweep U = U\ = U2 through J, 
comparing the RG with the alternative method outlined 
above. We see that the usual RG performs reasonably 
well, implying that the harmonic approximation is safe. 

Finally, we consider another potentially dangerous sce- 
nario. We return to the lattice shown in Figure 24. Now, 
we assume that J12 is greater than all other Josephson 
couplings, but it too is much weaker than the charging 
energies U\ and U2. In particular, jj^ = 0.05. Then, 




FIG. 25: In this test, J = J12 is assumed to be the strongest 
coupling in the system, but U = U\ — U2 may be of the 
same order. An effective coupling between sites 3 and 5 is 
calculated using two methods. One is the usual RG scheme 
used in this paper. Another is a hybrid exact diagonalization 
and RG scheme: The two-site problem of sites 1 and 2 is 
exactly diagonalized. Then, the resulting cluster is decimated 
away, and perturbation theory is used to calculate an effective 
coupling between sites 3 and 5. The two candidate values 
for the effective coupling J 35 are compared in the plot, as a 
function of jj . 



we sweep U\ such that it passes through a regime where 
\U\ — [/2I < J 12- The danger here is that the RG may ig- 
nore resonance effects associated with this region. Within 
the usual RG, sites 1 and 2 will be decimated in turn to 
give: 



^35. 



RG 



J: 



34, RG 



J13J12J25 
U1U2 



Ui 



(D2) 



(D3) 



where we ignore subleading corrections coming from po- 
tential applications of the sum rule, depending on the 
order of decimation of sites 1 and 2. To consider poten- 
tial resonance effects, we can also implement the same 
hybrid exact diagonalization and RG procedure that we 
used above. In Figure 26, we compare the two methods 
and find excellent agreement. 
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FIG. 26: In this test, J12 = 0.05 and U2 = 1. Hence, J12 is 
relatively quite weak. However, we vary £/i such that it passes 
through a regime where \Ui — L/2I < J12, where there may 
be a danger of resonance effects. We calculate the effective 
couplings J34 and J35 using two schemes. One is the usual 
RG scheme used in this paper. Another is a hybrid exact 
diagonalization and RG scheme: the two-site problem of sites 
1 and 2 is exactly diagonalized. Then, the resulting cluster is 
decimated away, and perturbation theory is used to calculate 
an effective coupling between sites 3 and 4 and between sites 3 
and 5. The effective couplings predicted by the two methods 
are compared in the plot, as a function of No resonance 
effects are observed. 



4. Strong Disorder RG vs. Exact Diagonalization 

As a final test of the RG procedure, we now compare 
the RG to exact diagonalization of small systems. We 
truncate the possible number fluctuation on each site to 
rij = —1,0,1, interpret these three values as possible z- 
axis spin projections of a spin-one object, and in doing 
so, arrive at a "spin-one" model: 



H 



(jk) 



Jjk 



{S+S^ + SjSD + Y.UjSl (D4) 



The Hilbert space of this spin-one model grows with the 
size of the lattice as 3 L . The particle number conser- 
vation of the rotor model manifests here as total spin 
conservation along the z-axis. This means that we can 
partition the Hilbert space into different total spin sec- 
tors and diagonalize the sectors separately. For most 
ground state expectation values, we just need to diago- 
nalize the total spin zero sector, and to calculate a charg- 
ing gap, the only additional diagonalization needed is for 
the total spin one sector. Despite these simplifications, 
computational limitations restrict us to studying 3x3 
lattices using CLAPACK. Testing the RG against exact 
diagonalization cannot directly tell us about the reliabil- 
ity of the RG at criticality, because exact diagonalization 
is limited to very small system sizes. However, a compar- 
ison with exact diagonalization can tell us how well the 
RG captures information about small patches of a larger 
lattice, and this information is potentially quite valuable 
for building confidence in the RG. 

Another complication arises precisely due to the 
Hilbert space truncation: the spin-one model may not 
always approximate the full rotor model well. This is 



FIG. 27: A sample- by-sample comparison of the particle num- 
ber variance predictions from exact diagonalization and from 
the renormalization procedure. The disordered couplings are 
drawn from three different choices of the distributions, with 
100 samples per distribution type. See the text of Appendix 
D for details on the distribution choices GG, GP, and FP. 
Also pictured is the coincidence line y = x, along which the 
points would ideally fall for full quantitative agreement. 
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FIG. 28: Same as Figure 27, but the quantity being compared 
is the charging gap A m i n . 



especially true when there is cluster formation, because 
then the rotor model has more of an opportunity to access 
particle number fluctuations that exceed 1 in magnitude. 
In other words, the strong disorder renormalization group 
and the exact diagonalization of the spin-one model con- 
stitute different approximations to the behavior of our 
random boson model. We cannot expect the two ap- 
proximations to show full quantitative agreement, but 
we proceed with the comparison, despite its limitations, 
with the hope of at least seeing qualitative correspon- 
dence between the two methods. 

Our general approach to comparing the RG with exact 
diagonalization will be to measure physical quantities, on 
individual 3x3 samples, using both methods and assess 
if there is a correlation between the predictions. The first 
quantity that we compare is the particle number variance 
(Bl). The interested reader may consult Appendix B to 
see how this quantity is calculated during the renormal- 
ization procedure. We also compare the charging gap 
Amin 5 the minimum energy needed to add a particle or 
hole to the system. This quantity is typically estimated 
during site decimation. The logic behind site decimation 
is that the charging energy for some site X is greater 
than all other scales in the problem, so the site can be 
disconnected from the rest of the lattice to leading or- 



31 



der. Then, the charging energy gives an estimate for the 
local charging gap on that site. During the RG, we find 
many such charging gaps from the various site decima- 
tions. The minimum among all of these gives an estimate 
for the charging gap for the whole system. This minimum 
is always given by the charging energy of the final remain- 
ing site, so an estimate of A m i n can be simply obtained 
by renormalizing down to a single site problem and mea- 
suring the charging energy of the remaining site. For the 
purposes of comparison to exact diagonalization however, 
we find that we obtain better quantitative agreement be- 
tween the two methods if we renormalize down to an 
effective two-site problem and then perform exact diag- 
onalization on that system. Exactly diagonalizing the 
n tot — an d ritot = 1 sectors of this two-site problem 
then yields a charging gap for the system. In this exact 
diagonalization, we need not truncate the on-site number 
fluctuation to tlx — —1,0,1. Instead, in the numerics, 
we typically truncate to nx = —100 . . . 100. 

In Figures 27 and 28, we present comparisons for three 
different data sets. In the first data set, we use take Pi(U) 
and Pi(J) to be Gaussian. We fix Uo = 10 and <jjj = 3. 
Then, we randomly sample Jo in the interval (0, 5) and 
a j in the interval (0, The aim is to approximate 
some of the environments that the RG encounters in runs 
such as those reported in Figure 6. The second data set 
uses the distributions described in Appendix C: P%(JJ) 



is Gaussian and Pi(J) oc J -1,6 . We randomly sample 
Uo G (6.5,20). Here, the motivation is to look at the 
types of environments that the RG encounters when it 
approaches the unstable fixed point from above. In the 
final data set, we try to mimic 3x3 patches that the RG 
might encounter near criticality. To this end, the initial 
J distribution is fixed to a power law P(J) oc J -116 
(see equation (39)) and the cutoffs are chosen so that the 
ratio of J m i n to J max is approximately that observed in 
panel (b) of Figure 13. The distribution P%(U) is flat with 
U max = Jmax and with f7 m i n randomly sampled such that 
the ratio of /7 m i n to U max lies in (e _2 ,e _1 ). Figures 27 
and 28 identify these three data sets with the labels GG, 
GP, and FP respectively. 

Both figures show that the predictions of the RG are 
clearly correlated with the predictions from exact diago- 
nalization, although the level of quantitative correspon- 
dence varies. For the particle number variance, quanti- 
tative agreement is lost at higher values of the variance, 
essentially corresponding to cases in which there is clus- 
tering. One potential source of error could be the Hilbert 
space truncation of the exact diagonalization, although 
the structure of the harmonic ground state (27) makes 
it unlikely that this could account for the entire discrep- 
ancy. Nevertheless, these comparisons suggest that the 
RG is retaining useful information about the system. 
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